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MODELING OF THREE-DIMENSIONAL MIXING AND REACTING DUCTED PLOWS 

By 

S. W, Zelazny, A. J. Baker and W. L. Rushmore 
Textron's Bell Aerospace 

SUMMARY 


A computer code, based upon a finite element solution algo- 
rithm, was developed to solve the governing equations for three- 
dimensional, reacting boundary region, and constant area ducted 
flow fields. Effective diffusion coefficients are employed to 
allow analyses of turbulent, transitional or laminar flows. The 
code was used to Investigate mixing and reacting hydrogen jets 
Injected from multiple orifices, transverse and parallel to a 
supersonic air stream. Computational results provide a three- 
dimensional description of velocity, temperature, and species- 
concentration fields downstream of Injection. 

Experimental data for eight cases covering different Injec- 
tion conditions and geometries were modeled using mixing length 
theory (MLT) . These results were used as a baseline for examining 
the relative merits of other mixing models. Calculations were 
made using a two-equation turbulence model (k+d) and comparisons 
were made between experiment and mixing length theory predictions. 
The k+d model shows only a slight Improvement in predictive ca- 
pability over MLT. Results of an examination of the effect of 
tensorlal transport coefficients on mass and momentum field dis- 
tribution are also presented. Solutions demonstrating the ability 
of the code to model ducted flows and parallel strut Injection 
are presented and discussed. 


INTRODUCTION 


The hydrogen- fueled scramjet engine Is a prominent candidate 
for propulsion of advanced hypersonic cruise vehicles. (See for 
example, Becker and Klrkham (ref. 1) and Bushnell (ref. 2).) An 
airframe-integrated underbody engine configuration (figs. 1(a) 
and (b)) has been suggested (ref. 3), and design considerations 
are discussed by Henry and Anderson (ref. 4). Many alternative 
scramjet designs have been proposed by the U.S. Air Force, the 
U.S. Navy, and NASA. In all cases, however, fuel Introduction 
typically consists of rows of circular, choked-flow fuel Injector 
orifices mounted flush and normal to the combustor wall or In fins 
spanning the combustor Inlet (see fig. 1(c).) The various proposed 


component designs have largely emerged from laboratory experimen- 
tation wherein empirical relations have established a preliminary 
configuration. Detailed experimental parametric evaluations are 
then utilized to optimize design configuration. 

The ability to analytically predict turbulent, mixing, and 
reacting three-dimensional flows, and hence avoid the more costly 
exclusively experimental approach, has been the long-range goal 
of rocket and ramjet designers for more than a decade. Three 
very difficult problems must yield to solution to attain this 
goal. First, a computational technique for solving the appropri- 
ate three-dimensional flow field with a predominant flow direction 
is required. Second, proper turbulent diffusion models must be 
selected or developed, since the accuracy of the predictive cal- 
culation is strictly dependent upon the adequacy of these models 
for combined laminar and turbulent diffusion of mass, momentum, 
and energy. Finally, detailed baseline data characterizing the 
flow phenomena over a reasonably wide range of flow parameters 
must be obtained to confirm the validity of flow theoretical mod- 
eling . 

The objectives of this investigation were to - (1) use the 
COMOC code developed in earlier studies (refs. 5 , 6, and 7) to 
investigate the accuracy of existing mixing models in character- 
izing scramjet combustor flow fields, (2) develop a new mixing 
model if required, and (3) extend the applicability of the COMOC 
code to consider constant area ducted flow fields. Although this 
investigation deals primarily with scramjet oriented problems, it 
is Important to note that the analysis may be employed in other 
problem areas characterized by three-dimensional reacting flows, 
e.g., chemical lasers, smokestack emission and gas turbine 
combustion chambers. 

Governing conservation equations, effective diffusion coef- 
ficient models, and methods used to compute reacting flows and 
pressure variations along the duct are presented in the METHOD OF 
ANALYSIS section. Details on the structure and methodology of 
COMOC are then presented in FINITE ELEMENT SOLUTION ALGORITHM. 
Current combustor design concepts (ref. 4) for scramjet engines 
employ fuel Injection both from transverse wall Injectors and from 
Internal struts containing both parallel and transverse injection 
orifices. Experimental data characterizing these types of condi- 
tions have been reported in references 8, 9 , 10, 11 and were 

modeled using the COMOC computer code. A brief description of 
these experimental studies and the details of the analytical 
modeling of this data using COMOC are presented in the NUMERICAL 
RESULTS section. Key results and recommendations for future ef- 
forts are presented in CONCLUDING REMARKS. Finally, an Appendix 
is Included where the required input and output control for the 


2 



COMOC code Is described in a card by card sequence. 

NOMENCLATURE 


a 

A 

b 

B 


c 

c 


P 


c 


e 


C 

C 

G 

d 

f 

g 

h 

H 

i 


f 

k 


i j j ,k 


k 

K 

i 


boundary-condition coefficient 

species; area; Van Driest Damping factor 

coefficient 

species 

coefficient 

specific heat 

ratio of 

species 
skin friction 

mixing length constant, Eq . (13) 

differential; orifice diameter; turbulence dissipation, 
Eq. (16) 

function of known argument 
function of known argument 
static enthalpy; duct height 
stagnation enthalpy; hydrogen 
index 

unit vectors of rectangular Cartesian coordinate system 
thermal conductivity; turbulence kinetic energy 
generalized diffusion coefficient; equilibrium constant 
differential operator; number; mixing length 
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L 

m 

M 

n 


N 

N 

N 

N 

N 

N 

N 

0 

P 

q 


d 

k 

Le 

Pr 

Re 

Sc 


a 


r 


Q 


R 


s 


S 


T 


u . 

1 


u,U 


— 

u. u . 

1 J 


characteristic length; differential operator 
number 

Mach number; number of finite elements 

unit normal vector; number; nodes per element; 
dimensionality 

nitrogen; composition matrix 
dissipation number = 

turbulence kinetic energy number = 

Lewis number 

Prandtl number 

Reynolds number 

Schmidt number 

oxygen 

pressure 

generalized dependent variable 
dynamic pressure ratio 

generalized discretized dependent variable 

domain of elliptic operator; universal gas constant 

injector spacing 

mass source term; finite element assembly operator 
temperature 

velocity component in 1th direction 
velocity 

turbulent shear stress component 
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friction velocity, / t ’/ p 

w 

molecular weight 

rectangular Cartesian coordinate system 

= UgX^/V 

species mole fraction 
species mass fraction 


species 1 
pressure 
ratio of 
closure o 
boundary- 
increment 
kinematic 
kinematic 

kinematic 


dentlflcatlon 

gradient parameter; elemental species 
specific heats 

f elllptlcally coupled solution domain 
layer thickness 

eddy viscosity 

eddy viscosity resulting from u^u' 
eddy viscosity resulting from u^u^ 


turbulence dissipation diffusion coefficient 


turbulence kinetic energy diffusion coefficient 


mixing efficiency 
coefficient 

multiplier; turbulence sublayer constant 

viscosity 

density 

Integral kernel 



T 


integral kernel; wall shear 




equivalence ratio. 


(pu-, ) 

jet 

0.03^'Cpu^) 

■air 


$ functional 

X domain of initial-value operator 

w turbulence damping factor 

global solution domain 
Superscripts : 

e effective value 

T matrix transpose 

a species Identification 

unit vector 

* approximate solution 

Subscripts : 


e 

i,j 

m 

o 

t 

T 


global reference condition 
local reference condition 
tensor indices 
mth subdomain 
initial 

stagnation or total 
turbulent 


Notation : 

{ } column matrix 

[ ] square matrix 
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u 

union 

n 

Intersection 

I 

summation 

e 

belongs to 


METHOD OP ANALYSIS 


Many researchers are now giving attention to numerical solu- 
tion of three-dimensional parabolic and/or boundary-region flow 
fields. Most procedures employ a finite difference solution al- 
gorithm for variously combined forms of the continuity, momentum 
and energy equations. Note that the three-dimensional boundary- 
layer equations result from this parabolic set for flow fields 
wherein diffusion In one direction only Is Important and the cor- 
responding pressure gradient Is negligible. Several researchers 
have obtained solutions for the three-dimensional boundary-region 
flow of single-species fluids. Pal and Rubin (ref. 12) employ 
asymptotic expansions of the flow variables for laminar Incompres- 
sible flow after transformation to modified streamfunctlon and 
vortlclty. Results of extending the theory to a compressible 
perfect fluid In physical variables are reported by Crescl et al . 
(ref. 13 ), who used an extension of the numerical technique common 
to boundary-layer solutions. Extension to handle streamwlse pres- 
sure gradients and refinement of the overall method are reported 
by Rubin and Lin (ref. l4). Caretto et al . (ref. 15) present a 
finite difference algorithm for solution of three-dimensional 
boundary-region flows with extension to the "parabolic’’ Navler- 
Stokes equations. The results of computations for transitional 
Internal flows In rectangular ducts are presented by Curr et al . 
(ref. 16 ). Refinement of the overall procedure with particular 
attention to solution of the parabolic Navler-Stokes equations Is 
given by Patankar and Spalding (ref. 17). The key feature of 
their theory is a procedure for splitting the pressure field com- 
putation such that a two-dimensional boundary-value problem results 
for pressure In the transverse plane coupled to an assumed uniform 
streamwlse pressure gradient computed from global continuity. The 
latter step is similar to methods employed for computations In two- 
dimensional hydrodynamics (ref. I 8 ). 

Characterization of an n-specles, three-dimensional boundary- 
region or parabolic flow field requires solution of (n-1) species- 
continuity equations In addition to those previously mentioned. 
Caretto (ref. 15 ) and Patankar and Spalding (ref. 17) Include 
results of a finite-difference solution of heat, mass, and momentum 
transfer In three-dimensional parabolic flows. Baker (ref. 19) 
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presents a finite element solution algorithm for multiple-species 
diffusion in supersonic, three-dimensional boundary-region flow. 
However, no general three-dimensional solution algorithm has been 
published which considers mixing and reacting three-dimensional 
confined flows . 

The system of partial differential equations governing such 
three-dimensional, confined unidirectional flows of a compressible, 
reacting fluid Is obtained as an approximation to the full three- 
dimensional Navler-Stokes equations. This approximation, referred 
to as the "parabolic Navler-Stokes equations," describes steady, 
confined three-dimensional flows wherein (1) a predominant flow 
direction Is uniformly discernible; (2) diffusion processes In 
the predominant flow direction are negligible compared with con- 
vection; and (3) no disturbances are propagated upstream, e.g,, 
recirculation is not considered. 

Presented below and In FINITE-ELEMENT SOLUTION ALGORITHM Is 
a description of the governing equations which are solved In the 
COMOC computer code. A number of key points should be noted: 

(1) Solution of the three-dimensional parabolic Navler- 
Stokes equations Is obtained wherein the streamwlse pressure var- 
iation Is computed using the approach presented In the subsection 
Pressure Variations In Ducted Plows. This version of the com- 
puter code is referred to as C0M0C-3DPNS. 

(2) If a pressure distribution is known a priori then that 
section of the code which computes streamwlse pressure variations 
may be bypassed. The resulting equations are herein referred to 
as the three-dimensional boundary region equations and represent 

a subset of the 3DPNS system. This variant of the code Is referred 
to as C0M0C-3DBR. 

(3) The solution algorithm embodied In both the 3DBR and 
3DPNS codes Is limited to constant area flow domains. To con- 
sider variable area flows would require either a coordinate trans- 
formation, e.g., streamfunction or a method which would add finite 
elements to the -outer edges of the computational domain. Including 
this feature Into COMOC was considered beyond the scope of this 
investigation . 

The velocity vector lying on a three-dimensional Euclidean 
space spanned by a rectangular Cartesian coordinate system xq Is 
identified as 


/N /N /N 

Ui = u^l + U2J + u^k 


( 1 ) 
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For development of the governing equation system, assume that 1 
Is parallel to the predominant flow direction. Identify the two- 
dimensional vector differential operator 

( = J ( )>2 + ^>3 ^ 2 ) 

where the comma identifies the gradient operator. In Cartesian 
tensor notation with summation over 2 and 3 for repeated Latin 
subscripts, the parabolic Navler-Stokes equation system for a 
multiple-species, compressible, reacting flow takes the form 


0 = 


(puj_), + (pu^). 


(3) 




— !i Y? 

“Ic^Re 


pu^Y,k + S 


(4) 




,, u . I - pu, u . - P J . 


(5) 


k 


pu^H, . 




^^Re^Pr 


H, 


k),^ P""k^’k - 


1-N„ e 
Pr y 


Npr 




y rC 6 a j 

NscNpr %e ^ 


- N„ ,,e 
Sc Pr y 


( 6 ) 


The variables appearing in equations (3) to (6) are non-dlmen- 
sionallzed with respect to y^, p^, U^, H^, and a length constant 
L, and have their usual Interpretation in fluid mechanics. The 
Reynolds number , effective Prandtl number Np^, and effective 

Schmidt number N® are defined for a combination of laminar and 

Sc 

turbulent contributions as. 
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(7a) 


N 


Re 




L 

00 


TJ 


CO 


y 

N 

y 

N 



+ 


pe 


"Pr ' 


T 


pe 


^Sc 


(7 b) 
(7c) 


In equation (7)j y is the laminar viscosity, e is the kinematic 
eddy viscosity, and the subscript T denotes a turbulent reference 
parameter. The stagnation enthalpy Is defined In terms of species 
static enthalpies as 


H E ^ h“Y^ + ^ \ (8) 

a 

The static enthalpy Includes the heat of formation of the 
species In its definition as 


h 


a 


dT + h“ 

T 

o 


(9) 


An equation of state Is required to close the system. Assuming 
perfect-gas behavior for each species, from Dalton's law. 


p = PRT I ^ (10) 

a V/ 

vjhere R Is the universal gas constant and Is the molecular 
weight of the ath species. 

Turbulence Models 

Closure of the governing equations requires Introducing 
relationships to define effective viscosity, and turbulent Prandtl 
and Schmidt numbers. The momentum diffusion (shear stress) In 
the x^_ direction may be expressed as 


10 


y + PEi 2 )^ 1 , 2 ^) ,2 P^13^^1,3^) 


3 


( 11 ) 



where the boundary layer approximation has been used to Justify 
dropping (y + p£,,)u, -, ) . Three different approaches to 

-L± -L j -L j ^ 

modeling and were considered, (1) simple mixing length 

theory (MLT) where It was assumed e ^2 “ = Ej (2) MLT where 

£^2 ~ £ E ]_3 “ where Is a constant, and (3) the two- 

equation turbulence model reported by Launder and Spalding (ref. 
20). For MLT the kinematic eddy viscosity Is given by 


The mixing length Is defined as 


C, x„o) 
k 2 


X6(jj 


where 


= 0.435 

X = 0.09 

6 = boundary-layer thickness 

X 2 = coordinate normal to wall 


The Van Driest damping coefficient Is 



where 


» u»X2 

^2 = “v 

* . 

u = friction velocity, /t/p 

T = skin friction 

p = density at wall 

V = kinematic viscosity 

A E 23.5 


(0 < X 2 < A6/C^) 

(X 2 > X6/C^) (13) 
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The second type of model considered wherein e ^2 ^ 

Interest since the effects of cross flow varlati ons In on flow 
field development may be examined- Of particular Interest Is the 
case where ^ e ^2 since here the effect of modeling each com- 

ponent of the turbulent shear stress as a tensorlal quantity may 
be examined. 

.The final model considered, l.e., a two-equation turbulence 
model, ref. 20, p. 97, requires defining the scalar kinematic 
viscosity, e, In terms of a characteristic velocity and length 
scales as 

e = (1^1) 


where the velocity scale Is related to the turbulence kinetic 
energy k given by 


k 





(15) 


and the turbulent length scale, Is related algebraically to 

the turbulent dissipation, d, through equation (l6). 




k3/2/d 


( 16 ) 


The turbulent kinetic energy and dissipation are computed 
from solution of equations (17) and (l8) 


pu^k,^ 


pe 


N, N 




k Re 


pu^k.i + ££ 


Re 


(u^^2^ - ^17) 


pu-i d. 


l^’l 


PE ^ 

I d Re ! 


C^d 

Re 


(ui 2 )^ - C 2 P ^ (18) 


where It was assumed that the production of turbulent kinetic 
energy Is due to the u^u^ component of shear stress, l.e., u^_ ^ 

Is negligible with respect to u^ 2 ' This assumption Is reasonable 
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for the flow regions examined herein using the k and d model. 

The constants used In equations (1^), ( 17 ), and (l8) as suggested 

In ref. 20 are as follows: = 0.09, = 1.44, = 1-92, 

= 1.0, and 'N^ = 1.3. 


In addition to having to Introduce an expresslon(s ) to com- 
pute the effectljre kinematic viscosity, it is also necessary to 
define the turbulent Prandtl and Schmidt numbers. In this in- 
vestigation it was assumed that mass and thermal energy, both 
scalar quantities, diffuse at the same rate, hence the turbulent 


Lewis number, (N, ) , is equal to unity. The validity of this 

1j0 rji 

assumption is supported by numerous experimental studies and is 
generally used in most turbulence modeling efforts, e.g., see 
ref. 21. This assumption requires (Np ) = (No ) since (Ny ) = 


complete 

Generally 


, hence the only remaining information required to 
T 

the description of the turbulence parameters is (Np ) . 

^ T 

, it is sufficient to define (Np ) equal to a constant 

T 


in the range from 0.7 to 1.0. The results shown in NUMERICAL 
RESULTS were obtained with (Np ) = 0.9, 0.7 or using an empirical 

T 


correlation . 


Pressure Variations in Ducted Flows 


For internal flows, characterized by boundary-layer thick- 
nesses which are small in comparison with the overall internal 
duct dimension, the pressure distribution can be accurately ap- 
proximated by invlscld flow solutions. However, in the alternate 
case where the flow is confined in a duct whose lateral dimension 
is not large with respect to the boundary-layer thickness, this 
approach is invalid. Here, boundary-layer development directly 
influences the pressure distribution within the duct, and an axial 
pressure gradient is induced by viscous effects. For these flows; 
a quasl-one-dlmenslonal integral treatment of equations (3) and 
(5) has been suggested (refs. 17 and 21), wherein for steady 
flows, equations (3) and (5) are integrated across the duct trans- 
verse dimensions to obtain an equivalent expression written on 
mass-averaged dependent variables defined by 


Q 


1 

mA 


A(x^) 


pu^qdT 


(19) 
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In equation (19) j A(x ) is the duct area, which may be a function 
of axial location x^; iTi Is the mass flow rate 

rti = pu^A (20) 


and q represents a generalized dependent variable which may be 
selectively streamwlse velocity u^, static temperature T, or 

density p. Taking the logarithmic differential of equation (20) 
gives 

dA _ dit '^’^1 dp 

A m u^ p ^ 

The Integral momentum equation (eq. (5)) Implies 


Adp + Fdx^ + u^diti + itidu^ = 0 


( 22 ) 


where F Is the retarding force per unit length of duct exerted by 
viscous Interaction of the confined flow with the wall. The equa- 
tion of state for a perfect fluid of constant molecular weight may 
be logarithmically differentiated to yield 


iP = ^ ^ 

P p T 


Combining equations (19) to (23) yields an explicit relation for 
axial pressure gradient as 


P 


1 


F 

A 




mu^ 


Ap 


mu-, 

— ^ T 
AT ’1 


(24) 


* 

Area must be constant In any application of C0M0C-3DPNS since 
the solution algorithm for the governing equations (3,4,5 and 
6) Is limited to this condition. However, the solution algo- 
rithm for streamwlse pressure gradient has been formulated such 
that variable areas may be considered once the constraint on 
Eqs. (3), (4), (5), and (6) Is removed. 
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If an initial pressure level and the detailed flow field at 

a given station are known, equation (2H) can be evaluated and 

Integrated to yield downstream pressure levels. To achieve this, 

the friction force per unit duct length is related to the wall 

shear stress t as 
w 


F 


T dL 
J w 

s 


(25) 


where the line integration is performed about the boundary of the 

computational domain, e.g., about symmetry planes t = 0, whereas 

w 

along walls x is finite. For flow fields which may be charac- 
terized by two parallel walls separated by distance h and two 
symmetry planes separated by distance h, the frictional force per 
unit area would be given by 


F 

A 


'U 


T dL 
w 


X dL 
w 


/(Lh) 


(26) 


where the integrations are performed over the upper (U) and lower 
(L) walls and where x^ is evaluated as a function of the local 

flow properties near the wall (ref. 21) from the expression 


X = C, ^pu? (R*^ - 0.156R?'^^ + 0.08723 r2'^ 

w k 1 " * * 

+ 0.0371R;°’^^)F (27a) 

* P 

where C is the same constant (von Karman ' s ) used in defining the 

^ P 'V ^ 

mixing length in MLT, R^ = k R and R = 11 ^X 2 /^ • symbol (~) 

refers to parameters evaluated in the constant shear stress region 
hence the length scale, X 2 represents the distance above the wall 

where convective flow effects are negligible. Similarly, u^ 
represents the x, velocity component at the edge of the constant 
shear stress region. The effect of pressure gradient on shear 
stress is Included in equation (27a) by the parameter F given by 


1 - 




7i2 . 8^ ' ^ R J' ^)^ • 


1.6 



(27b) 
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Calculation of the change in mass flow rate with respect to 
axial distance requires a computational distinction between the 

3? f* 

actual mass flow it and the computed mass flow m . The difference 

. r . f 

between m and m provides an estimate of the pressure and pres- 
sure gradient to maintain conservation of mass. The rate of change 
of mass flow vrlth respect to Is defined as 


II 

1 — 1 

r\ 

Am 

Ax^ 

(28) 

where 



Am = 

f 

Am"^ - Am"^ 

( 29 ) 

and 



Aih^ = 

A^(x^ + Ax^) - m^(x^) 

(30a) 

II 

•s 

<3 

A^(x^) - m^(x^) 

(30b) 

II 

i — 1 

X 

<1 

(x^) - (x^) 

new old 

(31) 

In equation 

(30a), Am^ represents the mass 

flow Increment which 


f 

results from the mass addition, whereas Am (equation (30b)) 
represents the mass flow error obtained at the upstream station 
x^. Examination of equations (24) and (28) to (30) shows that 

use of equation (28) will always provide a pressure gradient which 
tends to make the computed and actual mass flow discrepancy decrease 
and hence model the physical flow. 

Reacting Flows 

There are two approximate methods which may be effectively 
used to describe reacting hydrogen-oxygen-alr systems. In the 
first case, assume that prototype scramjet combustors are ade- 
quately described by equilibrium combustion. The following 
reactions are operative: 

2H + 0 i H2O 

2H i 
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20 t O 2 
H + 0 i OH 

N 2 + 20 t 2NO ( 32 ) 


The equilibrium composition of the combustion byproducts is de- 
termined by applying the law of mass action to each reaction 
defined In equations (32). This yields definition of a set of 
equilibrium reaction constants K, which, for the simple reaction 

-y 

nA + mB C , are expressed In terms of species mole fraction X 

as 


K 


r.^Ai 

in 

r^Bi 

X 

1 



m 


(33) 


Solution of equation (32), coupled with equation (33) and con- 
servation of total and elemental mass, yields after linearization 
algebraic equation system for determination of the equilibrium 
composition of the system, of the form. 



{x“} 


{Constant} 


(34) 


In equation (34), the elements of the matrix [n] account for the 

ct 

distribution of the particular species mole fraction {X } con- 
taining the 3th elemental material, for example, 0, H, and N. 
Solution of the equilibrium temperature and species concentration 
requires an Iterative solution to a nonlinear algebraic equation 
system. 

A considerably less expensive method (from the standpoint of 
computer time) may be employed to obtain an upper limit on the 
effects of heat release on the flow field development. The mole 
fractions of the dissociated species 0 and H are usually small 
compared with those of 0 ^ and H^. Equations (32) may then be 

considerably simplified by assuming that the complete reaction 


H 


2 


+ 



H^O 


(35) 


Is the only reaction that occurs. In this case all the reacts 
with the available O 2 to form H 2 O. By describing the variation 
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of specific heats with temperature through a relationship of the 

form c = a + bT, the temperature is solved explicitly in terms 
P 

of the enthalpy and pressure without iteration. 


FINITE ELEMENT SOLUTION ALGORITHM 


The parabolic Navier-Stokes equation system and the three- 
dimensional boundary-region equation system excepting global con- 
tinuity (equation (3)) are uniformly constituted as initial bound- 
ary-value problems of mathematical physics. Each of the subject 
partial differential equations (equations (4) to (6)) is a special 
case of the general, second-order, nonlinear partial differential 
equation . 


L(q) 


K(q)q, 


k 


+ f (q,q, . ,x ) + g(q,x) 

k , 


0 (36) 


where q is a generalized dependent variable identifiable with each 
computational dependent variable. In equation (36), f and g are 
specified functions of their arguments, x is identified with 

for parabolic flows, and x^ are the coordinates for which second- 

order derivatives exist in the lead term. The finite-element 
element solution algorithm is based upon the assumption that L(q) 
is uniformly parabolic within a bounded open domain that is, 
the lead term in equation ( 36 ) is uniformly elliptic within its 
domain R, with closure 9R, where 


^ = R X [Xq,x) (37) 

and Xq i Xj is the range of the initial value operator. 

Equation Development 

If equation (36) is uniformly parabolic, unique solutions for 
q are obtained upon specification of functional constraints on 
9q = 9R X [XqjX) an Initial-condition specification on 

RU9R X . For constraints on 9q, the general form relates the 

function and its normal derivative everywhere on the closure 9R 
as 


A(q) = a^^^q(x^,x) + a ^ ^ ^Kq ( x^ , X ) , ^^n^ - a^^^ = 0 (38) 
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In equation (38), the a^^^(x^,x) are user-specified coefficients, 
the superscript bar notation constrains to 8R, and n^^ is the 

local outward-pointing unit normal vector. For an initial dis- 
tribution, assume that 


q(xj^,Xo) = 


(39) 


is given throughout RU9R x Xq • 

The finite-element solution algorithm is established for 
the equation system (28) to (39) by using the method of weighted 
residuals (MVJR) formulated on a local basis. Since equation (36) 
is valid throughout U, it is valid within disjoint interior sub- 
domains described by (Xj.,x)eRjjj x [xq,x)j called finite elements, 

wherein UR^ = R. An approximate solution for q within R^ x [XqjX)j 
called q*(x^,x)j is formed by expansion into a series solution of 
the form 


'^■*‘'^1"^^ = {®(x^)}^{Q(x)}j^ 


(^ 0 ) 


In equation (^0), the functionals subsets of 

set that is complete on Rj^. The expansion coefficients 
represent the unknown x~dependent values of Q*(Xj_,x) 
locations interior to R^^ and on the closure 9Rjjj, called 
the finite-element discretization of R. 


a function 

Qic(x) 

specific 

nodes of 


To establish the values taken by the expansion coefficients 
in equation (40), requires that the local error in the approximate 
solution to both the differential equation L(q^) and the boundary- 

condition statement ^(q^)? f’oi’ 9Rj^n9R 0, be rendered orthogonal 

to the space of the approximation functions. By employing an 
algebraic multiplier X the resultant equation sets can be combined 
as 


S 

m 


R 


{$(x. ) }L(q*)dx - 
1 m 


m 


9R, 


{il>(x. ) }£(q*)da = {0} (4l) 

H9R 1 ^ 


m 
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where S Is the mapping function from the finite-element subspace 
m 

to the global domain R, commonly termed the assembly operator. 

The number of equations (^1) prior to assembly is identical with 
the number of node points of the finite element R^. 

Equation (4l) forms the basic operation of the finite element 
solution algorithm and of the COMOC computer program. The lead 
term can be rearranged, and X determined by means of a Green-Gauss 
theorem: 


R 






m»k 


m 


’k 


dT = 


K p 


9R. 


}Kq*,j^nj^da 


m 


- K 


R. 


{$(x. ) } , Kq* , dT 
1 k m k 


m 


(42) 


For SRPiSRj^ nonvanishing (equation (42)), the corresponding segment 

of the closed-surface integral vrlll cancel the boundary-condition 

contribution (equation (4l)) by identifying Xa with k of equa- 
tion (36). The contributions to the closed-surface Integral 
(equation (42)), where 8R^n9R = 0, can be made to vanish (ref. 6). 

When equations (38) to (42) are combined, the globally assembled 
finite-element solution algorithm for the representative partial 
differential equation system becomes 


S 

m 


- K + g*)dx 


R 


m 


R 


■m <=m" 


- K 


I - a<3)) 

V m ^m m ' 


da 


^ aRjjjOaR 


= {0} 


(43) 


The rank of the global equation system (equation (43)) is identical 
with the total number of node points on RU3R for which the de- 
pendent variable requires solution. Equation (43) is a first- 
order, ordinary differential system, and the matrix structure is 
sparse and banded. Solution of the ordinary differential equation 
system is obtained by using a predictor-corrector finite difference 
numerical integration algorithm (ref. 6). 
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A solution algorithm Is required for the continuity equation, 
which Is retained as equation (3) for boundary-region flows. 

Since equation (3) Is an Initial-value problem on pu 2 as a func- 
tion of X 2 , with and x^ appearing as parameters, the approxi- 
mation function need span only the transverse coordinate direc- 
tion as 


(pup)* = {3>(xp) }'^{pV(x. ,Xp) } (44) 

The matrix elements of pV are nodal values of pu^; their functional 

dependence requires solution of equation (3) along lines (x^,x^) 

equal a constant. Since equation (3) exists In standard form as 
an ordinary differential equation, direct numerical quadrature 
yields the required solution at node points of the discretization. 

Structure of the COMOC Code 

The COMOC computer program system Is being developed to 
transmit the rapid theoretical progress In finite element solution 
methodology Into a viable numerical solution capability. In the 
course of generating this general-purpose system, several variants 
of COMOC have been developed for specific problem classes. In- 
cluding transient thermal analysis and the two-dimensional Navler- 
Stokes equations as well as the three-dimensional boundary-region 
equations. The present operational variant of COMOC Is capable 
of solving each of these problem classes and has been extended to 
Include the parabolic Navler-Stokes equation system. An on-line 
restart feature allows the user to switch between boundary-region 
and parabolic Navler-Stokes systems according to the requirements 
of the problem at hand. 

The finite element solution algorithm Is utilized to cast 
the original Initial-valued, elliptic boundary-value problems Into 
large-order systems of purely Initial-value problems. The program 
then Integrates the discretized equivalent of the governing equa- 
tion system In the direction parallel to the predominant flow. 
Initial distributions of all dependent variables may be arbitrarily 
specified, and boundary constraints for each can be specified by 
the user on arbitrarily disjoint segments of the solution domain 
closure. The solutions for each dependent variable, and all com- 
puted parameters, are established at node points lying on a 
speclflably nonregular, computational lattice formed by plane 
trlangulatlon of the elliptic portion of the solution domain. 

Each of the computational triangles Is spanned by a linear approx- 
imation function used for all Independent and dependent variables 
as well as each solution parameter. 
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The COMOC system Is built upon the macrostructure Illustrated 
in figure 3- The main executive routine allocates core by means 
of a variable dimensioning scheme based upon the total degrees of 
freedom of the global problem. The size of the largest problems 
that can be solved Is thus limited by only the available core of 
the computer In use. The precise mix between the number of de- 
pendent variables (and parameters) and fineness of the discreti- 
zation is user-specifiable and widely variable. The input module 
serves Its standard function for all arrays of dependent vari- 
ables, parameters, and geometric coordinates. The discretization 
module forms the finite element discretization of the elliptic 
solution domain and evaluates all required finite-element non- 
standard matrices and standard-matrix multipliers. The Initial- 
ization module computes the remaining initial parametric data 
required to start the solution. The Integration module consti- 
tutes the primary execution sequence of problem solution. It 
utilizes a highly stable, predictor-corrector integration algo- 
rithm for the column vector of unknowns of the solution. Calls 
to auxiliary routines for parameter evaluation (viscosity, Prandtl 
number, source terms, combustion parameters, etc.) as specified 
functions of dependent and/or independent variables are governed 
by the Integration module. The user has considerable latitude to 
adapt COMOC to the specifics of his particular problem by directly 
Inserting readily written subroutines to compute special forms of 
these parameters. The output module Is similarly addressed from 
the Integration sequence and serves Its standard function via a 
highly automated array display algorithm. COMOC can execute dis- 
tinct problems In sequence and contains an automatic restart 
capability to continue solutions. 

Accuracy and Convergence 

The three-dimensional boundary region equation system solved 
using COMOC may be routinely employed to consider two-dimensional 
problems as a special case. This feature Is Important since the 
COMOC generated results may then be evaluated for accuracy and 
convergence by comparison with solutions produced by finite-dif- 
ference techniques and with a similarity solution for constant 
specific heat. With this point In mind, three check cases were 
considered and are presented below. 

Compressible Boundary Layer - Consider a nominal Mach 5, laminar, 
two-dimensional, air boundary-layer flovf over an adiabatic v;all 
In a favorable pressure gradient. With the assumption of constant 
specific heat, the flow Is Isoenergetlc and It Is necessary only 
to solve the x^ momentum equation and the continuity equation. 

The initial distribution for longitudinal velocity u^ Is estab- 
lished from the similar solution for B = 0.5 and S = 0 of reference 
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24. The Initial distribution for U 2 is obtained iteratively, 
and Sutherland's law is employed to compute viscosity. 

The test case is initialized at = 0.03 m downstream from 

the leading edge. The boundary-layer thickness at this station, 

<S , is 0.0039 m, the local Mach number, M , is 3- 77, the Reynolds 
o ® 5 

number per unit length, is O .83 x 10 per meter, and the 

adiabatic wall temperature, T , is 1000°K. Shown in figure 4 are 
the COMOC computed distributions of skin friction, local free- 
stream Mach number, and boundary-layer thickness for the case of 
constant specific heat. These were obtained with two uniform 
finite-element discretizations corresponding to four and eight 
elements spanning the initial boundary-layer thickness. The in- 
put static pressure distribution Pg(x^) is also presented for 

reference. Only small differences, on the order of about 2 per- 
cent, exist between the two solutions, the finer discretization 
producing a slightly larger skin friction and smaller local Mach 
number. Superimposed in figure 4 for comparison purposes are 
the results for the similar solution (ref. 26 ) and a 20-zone 
finite-difference solution obtained with the Von Mlses coordinate 
transformation. Agreement among the four solutions is excellent 
(within 2 percent) for skin friction. The similar solution for 
M lies between the COMOC and finite difference solutions, and 

overall agreement is within ±3 percent 

Shown in figure 5 are computed velocity profiles at = 

22 . 7 , which is about midway through the presented solution domain. 
Shown for reference is the initial longitudinal velocity profile 
with the node locations of the four-element discretization super- 
imposed. Both COMOC solutions produce u^_ distributions that are 

slightly more concave upward in the midregion in comparison with 
the similar or finite difference solution. The eight-element 
COMOC solution lies closer to the similar solution in the region 
where the two finite element solutions differ. The finite dif- 
ference solution lies appreciably below both the COMOC and similar 
solutions near the freestream. The computed transverse velocities, 
which are also plotted in figure 5, show only slight differences 
between the two discretization solutions. The trends of the COMOC 
solutions are in excellent agreement with the established pro- 
cedures. This check case result establishes an accuracy assess- 
ment of the solution algorithm for the three-dimensional boundary- 
region equations. 

Developing and Developed Channel Flow - Other check cases used 
in the evaluation of the parabolic Navler-Stokes equation system 
have been examined for three different channel flow configurations. 
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Figure 6 summarizes the results for a nonreacting subsonic flow 
to evaluate the ability of the pressure solution algorithm (equa- 
tion (24)) to compute a constant streamwlse gradient. For the 
fully developed channel flow, streamwlse velocity and the pressure 
gradient are computationally maintained to within ± 2.5 percent 
of their initial values. The computations for developing channel 
flow correctly predicted the downstream distance required to attain 
fully developed flow; that is, COMOC predicted that the flow was 
fully developed at x,/h = 33 compared with x-, /h = 30 reported in 
ref. 25. 

Mixing and Reacting Channel Flow - An evaluation similar to that 
described above was performed to assess channel flow computations 
with heat addition. Conditions were selected such that in the 
initial portion of the flow, reaction of hot air v;lth cold hydro- 
gen Induces a favorable pressure gradient (heat addition in sub- 
sonic flow). Hovrever, after the available oxygen supply is ex- 
hausted, the continued mixing of the cold hydrogen with the heated 
combustion produces an overall temperature drop, and hence, an 
adverse pressure gradient. The computational results are sum- 
marized in figure 7; the trends are observed to have been correctly 
predicted by COMOC while maintaining conservation of mass to within 
± 1.0 percent, l.e., the correct increase and then decrease in 
velocity as the gas initially heats and then cools. 


NUMERICAL RESULTS 

The objective of this study is to provide a theoretical 
model describing the three-dimensional mixing and reacting flow 
fields v/hich characterize scramjet combustors. Theoretical pre- 
dictions will prove useful in providing guidance in selecting 
fuel Injection/strut geometry and estimating combustor size and 
performance requirements. Before any theoretical model is used 
to provide design criteria, it is essential to test the ability 
of the model to determine its reliability of predictions, i.e., 
how close are the predictions to experimentally observed trends 
and are predictions generally above or below experiment. Ulth 
this point in mind the COMOC code v;as exercised to: (l) examine 

the sensitivity of predictions to a number of key parameters, 

(2) determine the accuracy of predictions for a different non- 
reacting flow condition, and (3) model reacting flows. 

Sensitivity Study 

The sensitivity of the predicted flow field solution to 
a number of key parameters was studied and results presented 
below. Specifically, those parameters examined were discreti- 
zation, eddy viscosity model, tensorlal character of eddy 
viscosity, turbulent Prandtl number, transverse velocity. The 
detailed experimental results of Rogers (refs. 8 and 9) for 
the configuration Illustrated in figure 8 provide the necessary 
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data base for* comparison of predictions. Initial conditions for 
the predictions were established from these data, and the down- 
stream station at x^/d = 30 was selected as the initialization 

station. The original raw data consist of a single vertical tra- 
verse and three lateral traverses on the transverse plane at sev- 
eral stations. The measured hydrogen mass fraction distribu- 
tions appear of Gaussian shape; however, the symmetry plane of 
the data was variously displaced from the geometric symmetry 
plane. Although the entire flow field could be computed numeri- 
cally, the strong appearance of a data symmetry plane suggested 
establishment of a corresponding computational solution domain. 
Therefore, a cubic spline Interpolation program was applied to 
the ravr data program to establish the x^/d location of the data 

symmetry plane via a minimization criteria on the wings of the 
Gausslan-type distributions. The spline package then interpolated 
the raw data for hydrogen mass fraction and Uj and output the 

evaluation of the Interpolation polynomials at node points of 
the finite-element discretization of the transverse plane. A 
representative case of the spline-computed distributions of hy- 
drogen mass fraction Is shown In figure 9 In comparison with the 
spread and context of the experimental data. 

Although plots of the form of figure 9 are geometrically 
aesthetic, the transition from the Initial distributions and 
significant detail on solution accuracy and trends are better 

H 

obtained by plotting concentration profiles (x^/d against Y ) 
along planes x^/d = constant at each longitudinal station for 

which data measurements exist, or simply examining variations 

along the centerplane In the normal direction and along the wall 
H 

(Y against x^/d at x^/d = 0 and fixed x^/d) . Both of these 

methods of comparing predictions with data were used and examples 
for case 1 of Table 1 are shown in figures 10 and 11 (the solid 
curves shown in figure 11 represent a best fit to the data; wall 
values were obtained by extrapolation). 

Turbulent mixing length theory has been very successful In 
modeling various types of flow fields as shown by Launder and 
Spalding, ref. 20. Therefore, before considering more complex 
turbulence models, the MLT given by equations (12) and (13) was 
used to establish a basis for comparison. As shown In figures 

10 and 11, MLT and the assumption of Np^ = 0.90, = 1.0 

results In predicted Hp mass fractions In good agreement with 
the data tov/ards the outer edge of the mixing boundaries. All 
dependent variable boundary conditions were taken to be zero 
gradient at both the freestream and wall except for the velocity 
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components which were specified to be zero along the wall. 
Predicted mass fractions near the wall are greater (H.H% 

against 3-2% at x^/d = 60, x^/d = 0) than the experimental values. 

As will be shown below, this trend carried over for most of the 
other flow conditions studied (Table 2). Possible explanations 
for this discrepancy may be: (l) inaccuracies in the data In the 
vicinity of the wall, (2) Inadequacies of the mixing model, and 
(3) specifying u^ = 0.0 due to lack of experimental data for 

this velocity component (see discussion of case 1-4 below). 


As with any multidimensional computations in compressible 
viscous fluid mechanics. It Is important to establish a quanti- 
tative accuracy assessment. For the cold-flow configuration 
studied and reported herein, an accuracy measure of the adequacy 
of the employed discretization Is possible by determining the 
conservation properties of the solution. For the cold mixing 
case, the species-continuity equation for hydrogen mass fraction 
can be written in explicit conservation form. Integrating this 
equation over a three-dimensional control volume and using Gauss' 
theorem (ref. 5) determines that the total hydrogen mass flow, 

.H 


pu^Y dA, would be rigorously conserved by an analytic 


that Is, 

solution."^ COMOC evaluates this parameter at each output station 
by using linear finite element approximation functionals for 
each variable and performing the Integrations analytically. 
(Thus, the order of the evaluation Is consistent with that of 
the solution of the partial differential equations.) 


For case 1-1 of Table 1, a monotonically Increasing loss of 
hydrogen mass flow v/lth Increasing distance dovrnstream was com- 
puted; at Xj/d = 120, the computed loss equaled 8.8 percent of 

the mass flow computed at station x^/d = 30. The 100-element 

standard nonuniform discretization was refined by a factor of 2 
In each coordinate direction to produce 400 finite elements 
spanning R (see fig. 12, diagonals omitted), and the computation 
was repeated on 30 <_ x^/d £ 60. Over this Interval, the coarse 

discretization yielded a computed 5 percent loss In hydrogen 
mass flow. The fine discretization produced a modest variation 
In computed hydrogen mass flow over the initial Interval, with a 
1 1/2-percent net loss computed by x^/d = 60. The resulting 

detailed differences In computed distributions of hydrogen mass 
fraction are shown in figure 13- Above the peak and outside the 
near wall region, differences are undlscernlble on the scale of 
the plots. Within the near wall region, the maximum difference 
in computed hydrogen levels Is less than 8-percent, which compares 
favorably with the 10 to 20 percent spread of the "best symmetry" 
data. The computational expense of these comparison solutions 
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differed by over an order of magnitude. On an IBM 360/65, using 
no out-of-core devices for either case, the CPU time of the 100- 
element solution was about 250 seconds; on the same Interval, the 
400-element solution required 2600 seconds. Two parts constitute 
this Increase (see ref. 6); a factor of about 4, due to the fact 
that the element DO loops In COMOC were 4 times longer, and a 
multiplicative factor of about 2, due to Increased solution stiff- 
ness resulting from the refined grid Itself. The ability of 
coarse finite-element discretizations, using low-order functionals, 
to preserve adequate engineering solution accuracy appears a dis- 
tinct feature of the algorithm. 

It Is of Interest to examine the sensitivity of predictions 
to variations In the empirical constants used In the mixing model. 
Case 1-3 of table 1 was used to examine the effect of decreasing 
the mixing length constant X from 0.09 to 0.07 while reducing the 
effective Prandtl number from 0.9 to 0.7. In the outer region 
of the boundary layer e X , hence decreasing X will slow the 
momentum mixing since e Is Initially decreased by 40^. Similarly, 
the mass diffusion coefficient Is also decreased by 22 percent. 

As shovrn In figure l4, the direct effect on mass diffusion by 

0 

these changes In X and Np^ are to slow the mixing rate of Hp such 

that the wall concentration Is lower than case 1-1 by 10 percent 
(4.,4^ vs 4.0^ at the wall and 4.5% vs 4.9% at the peak). 

The correlation between predictions and data for the test 
case shown In figures 10 and 11 Is good except In the near wall 
region. For good agreement In the centroldal region of the 
hydrogen Jet, It Is necessary that the maximum hydrogen concen- 
tration remain off the wall. Therefore, either from three-dimen- 
sional effects or a complex turbulence Interaction between the 
Jet and the wall, there may be mechanisms In play capable of re- 
sisting the unidirectional trend of maximum diffusion to the 
wall. It has been hypothesized, as a result of the Initial 
studies (ref. 5), that the existence of a mass flux transverse 
to the main flow direction and along the plate surface might 
account for the exnerlmentally measured centroldal peak. Such a 
transverse mass flux could be Initiated by the displacement effect 
of the sonic hydrogen Jet Issuing transverse to the main flow, 
since In such an Interaction problem, the Jet appears to the 
mainstream flow In many ways similar to an Impervious body. Con- 
sequently, Immediately downstream of the transverse Jet, there 
must exist an approximately spheroidal fixed recirculation region 
near the wall, and a transverse mass flux would be required to 
alleviate a localized low-pressure area Just downstream of this 
bubble . 

This hypothesis was computationally evaluated In test case 
1-4 by Imposition of a small negative, transverse velocity dis- 
tribution beneath the measured hydrogen concentration maximum 
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at x^/d = 30. All other conditions were the same as those used 

in case 1-1. The magnitude and vertical extent of the Imposed 
transverse velocity on the symmetry center plane (x^/d) = 0 are 

shown In figure 15a. Figure 15b Illustrates the lateral spread 
of the Imposed transverse velocity distribution along a specific 
(x^, X 2 ) plane, l.e., X 2 /d = 1.0. As shown (figures 15a and b), 

the Ug velocity component is less than one tenth its peak value for 

x^/d ^ 5o. 

The influence of the Imposed transverse velocity on the 
predicted distributions of hydrogen mass fraction at the two 
downstream, data stations Is shown In figure 16 . The selected u 
velocity distribution Is observed to not significantly alter 
the mass fraction distributions above the peak but does substan- 
tially promote the existence of a local off-plate maximum In the 
centroldal region at x^/d = 60 . However, by the time the last 

data station Is reached (x^/d = 120), the Imposed transverse 

velocity distribution has been essentially dissipated, and the 
computed distributions of hydrogen mass fraction in the centroldal 
region are noted to revert to the form of the maximum existing 
at the plate surface. It can be concluded, therefore, that trans- 
verse mass flow Is probably of Influence In the near region down- 
stream of the point of injection. However, there is as yet some 
undetermined mechanism for maintaining the off-axis peak In the 
mass fraction distribution at stations far downstream. This 
undoubtedly points to some deficiency In the turbulent mixing 
model for this configuration and/or the experimental wall values. 

As recently pointed out by Launder, Reece, and Rodl (ref. 28) 
while effective viscosity models have led to satisfactory pre- 
dictions In two-dimensional flows, their use In three-dimensional 
flows with more than one significant m.ean velocity gradient has 
achieved only moderate success. Considering a Reynolds-stress 
closure model for the problem of Interest herein would be pre- 
mature due to overall complexity of the flow field (compressible, 
variable molecular weight, not well defined Initial conditions 
and in particular lack of turbulence measurements). However, the 
tensorlal character of the turbulent shear stress may be examined 
using the COMOC code by characterizing the eddy viscosity as a 
tensorlal quantity. Shown In figure 17a and 17b are the predicted 
H 2 concentrations obtained by assuming that the eddy viscosity 

characterizing diffusion in one of the directions normal to the 
main flow coordinate Is twice that of the eddy viscosity In the 
other direction, l.e.. In figure 17a, e ^2 “ ^ 2e, 

whereas In figure 17b, ~ ^13 “ seen from the 

results shown In figure 17a or 17b, it would not be possible to 
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Improve the level of agreement between data and theory by a simple 
adjustment (change In a single constant) In the eddy viscosity 
since neither case gave any significant Improvement In the pre- 
dicted concentration distribution. These results suggest that 

direct modeling of each shear stress component may be required 
but further comparisons with other data are required before any 
definite conclusions may be drawn. 

The two-equation turbulence model (k+d) of Launder and 
Spalding (ref. 29) has enjoyed considerable success In modeling 
two-dimensional turbulent flows. In this study Initial values 
for k and d were obtained by computing an eddy viscosity from 
MLT, equation (12), mixing length from equation (13) » and using 
the relationship between k, d, 5., and e given by equations (l4) 
and (l6). As shown In figure l8, a modest Improvement In the 
correlation between theory and experiment was obtained using the 
k+d model of ref. 29 to predict the H 2 concentrations. To be 

noted, however. Is that the Introduction of two additional equa- 
tions Increased computer run time by approximately 30 percent 
and leads us to conclude that use of higher order turbulence 
models for the problem of Interest herein are somewhat premature 
at this time. 


Computations for cases 1-1 and 1-6 were made assuming a 
constant effective Prandtl number. However, a number of studies, 

e.g., ref. 30, have shown that , and hence, Np^ vary con- 

siderably across the boundary layer. Case 1-7 was considered to 
examine the effect of considering Prandtl number variations. 
Following Wassel and Catton (ref. 30), the expression given by 
equation (^5) for the turbulent Prandtl number was used to model 
conditions corresponding to q^ = 1.0, s/d = 12.5. 


N 


Pr) 


T 


1 - exp 
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1 Pr 
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Npr ( Pijt/ P ) 
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where C 2 _ = 0.21, C 2 = 5.25, C^ = 0.20, = 
values of viscosity and Prandtl number were 
l.e., p = 0.000018 N-sec/m^ and = O. 69 . 
Prandtl number Is obtained from equation 7b 


5 . 0 , and the laminar 
assumed constant. 

The effective 
and is given by 
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Predictions of the H 2 concentrations obtained using this model 

are shown in figure 19 where It Is noted that by considering a 
variable Prandtl number the rate of diffusion to the near 

wall region Is decreased in agreement with the experimentally 
observed trends. This result suggests that the optimum mixing 
model for the flow field of Interest would be MLT coupled with 
the variable Prandtl number expression given by equations (45) 
and ( 46 ) . 

Modeling Flows with Various Values of and s/d 

It Is essential that a number of different flow conditions 
be examined when evaluating any empirical model before drawing 
any conclusions about Its ability to model a class or classes of 
flows. Table 2 lists eight cases modeled using MLT and a constant 
effective Prandtl number; results are shown In figures 20 to 26. 

The performance of the k+d and variable N© models for these 

cases may be estimated by examining predictions obtained using 
these models for case 2-1 (figures 17 and I 8 ). In general, the 
level of agreement between predictions and experiment for these 
seven cases ( 2-2 to 2 - 8 ) are similar to those obtained In modeling 
case 2-1. Specifically, the predicted H 2 mass fractions were 

greater than the experimental values In five of the seven cases. 

To be noted is that mixing was predicted to be faster than observed 
for the three = 0.5 cases. A means for accounting for this 

effect through the mixing model would be to employ the k+d tur- 
bulence equations and relate the Initial turbulence kinetic 
energy to the Injection parameter, q , l.e., as q^ Increases 

the turbulence level In the boundary layer Increases. The com- 
puter code was not exercised to further explore this concept 
since the results would represent only a data fit and not provide 
any new Information. Detailed turbulence data would be required 
to test the ability of the k+d model to characterize the dominant 
turbulent mixing processes. However, from the comparisons between 
data and theory. It Is concluded that the COMOC code even with 

the simplest turbulence model considered (MLT and Np^ = 0.9) 

gives correlation with experimental data with sufficient accuracy 
to provide useful engineering design guidance. In the following 
subsections, we present results obtained using COMOC to characterize 
a number of scramjet Injection modes. 
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simulation of the Near Injection Region 

A distinct feature of the three-dimensional boundary-region 
and parabolic Navler-Stokes equation systems Is the capability 
to obtain a three-dimensional solution while marching In one 
coordinate direction. Elimination of the requirement for a down- 
stream boundary condition Is particularly Important. For the 
subject mixing and combustion studies, however, the corresponding 
penalty Is that an accurate Initial condition Is required to model 
adequately the complex near-lnj ectlon flow field under study. In 
the previous section, the computations took advantage of detailed 
experimental profiles for the establishment of Initial conditions. 
In the more general case, and In particular for hot-flow cases 
with combustion, detailed distributions of Initial conditions are 
specifically unavailable, and a theoretical device for establishing 
the starting point of the solutions Is required. Flow fields In- 
volving the parallel Injection of dissimilar fluids present no 
difficulty, since smooth transitions occur and boundary-layer and 
shear-layer concepts are appropriate. However, for transverse In- 
jection, this Is not the case, and some alternative means Is 
required . 

In anticipation of this need, a task In the early phases of 
the study was to evaluate the concept of a numerical "virtual 
source" as a means for eliminating the requirement for detailed 
Initial conditions (ref. 5). Injection of a jet from an orifice 
In a plate transverse to a supersonic alrstream has been the sub- 
ject of a number of Investigations. The Important correlating 
parameter appears to be dynamic pressure ratio q^. Most exper- 
imental data are for large values of q^, for which the Jet has 

sufficient momentum to penetrate the boundary layer and produce a 
complicated separation region and bow shock ahead of the jet. 
However, for the present cases, q ranges betv/een 0.5 and 1.5 3 

and a significant part of the jet remains within the turbulent 
boundary layer. Hence, mixing Is Initiated Immediately downstream 
of Injection. From these considerations, a theoretical model was 
proposed for establishing a barrel-shock - Mach disk hypothesis 
for turning of the transverse jet parallel to the main flow (see 
fig. 27 ). An analysis based on one-dlmenslonal considerations 
was developed to characterize the jet turning. The Important 
parameters In the model are dynamic pressure ratio q^ and free- 

stream Mach number and the output Is Injectant momentum and 

flow area. Details of the model are presented In reference 5- 

The validation of the concept of the virtual source as an 
Initial-condition generator was accomplished by using the detailed 
cold-flow experimental data discussed In the previous section. 
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The detailed predicted distribution of hydrogen mass fraction are 
shown In figures 28 and 29 for the virtual-source simulation of 
the standard test case (q^ =1.0 and s/d = 12.5) • For these 

0 

computations, transverse velocity u^ was assumed zero, Np^ = 0.9 

and MLT was employed. Even after marching only 30 diameters 
downstream from the point of Injection, agreement between the. 
predictions and the data (fig. 28) Is admirable, especially In 
the centroldal region, where there Is an excellent prediction of 
an off-plate peak. At the final station x^/d = 120 (fig. 29), 

agreement between the virtual-source simulation and data Is 
excellent, being essentially Identical with the results starting 
with data at x^^/d = 30 (fig. 10). 

Reacting Plows 

Three different reacting flow configurations were considered 
In this Investigation: (1) perpendicular Injection of Into a 

supersonic stream (either air or vitiated air) from a row of 
circular orifices aligned across a flat plate, (2) perpendicular 
Injection of from a circular orifice positioned on a strut In 

a scramjet combustor, and (3) tangential Injection of Hp also 

positioned on a strut In a scramjet combustor. The first case 
considered, referred to as case 3-1 In table 3, corresponds to 
case 2-1 using the virtual source simulation with the exception 
that here the flow Is allowed to react. The conditions for this 
reacting flow case are not significantly different from those 
reported by Rogers and Eggers (ref. 10) who have experimentally 
Investigated the reaction of hydrogen In a hot supersonic test 
gas flowing through a two-dimensional duct. (Data point ^ of 
reference 10 corresponds most closely to the conditions of case 
3-1.) Case 3-1 has been Introduced to examine the effects of 
heat release on the predicted flow field. Case 3-2 Is Identical 
to case 3-1 except the air stream has been replaced by a vitiated 
air stream. 

Reacting flow data for the strut Injection geometries have 
been reported by Anderson and Gooderum (ref. 11). These data 
(two cases) were simulated using the COMOC code, cases 3-3 and 
3-4 of table 3* The perpendicular Injection condition was simu- 
lated using the virtual source concept and the effects of shocks 
were neglected In each case except through Imposition of the 
reported static pressure variations with distance downstream. 
Cases 3-5 to 3-7 examine the effect of modeling reacting flows 
using the k+d model, and. In particular, the effect of turbulence 
on reaction. Cases 3-8 and 3-9 consider the effect of ducting 
the flow on the mixing. 
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H 2 Injection from a Plat Plate - Shown In figure 30 Is an alter- 


native method for presenting these data which was extended to 
reacting flow comparisons. This figure presents the peak hydro- 
gen concentration as a function of distance downstream. The 
agreement with data Is excellent in the range 30 x^/d <_ 120. 

The disagreement at x^/d = 7 is not serious In light of the 
further downstream agreement; the Indicated data point may well 
be significantly In error. Shown also In figure 30 are the tra- 
jectories of the peak hydrogen concentration above the plate 
x^/d and the lateral spread of the jet determined at the x^/d 

coordinate where the local hydrogen concentration equals 10 
percent of the local maximum. Observe that the local peak of 
the elevation trajectory sinks to the plate surface downstream 
of x^/d = 60 , as was observed for the eight cold-flow data cases. 

The plotted spread of computations In this region Indicates the 
span of X 2 /d for which the predicted hydrogen concentration varies 

by 10 percent. The agreement of lateral spreading rate with data 
Is excellent. 

At the lower right In figure 30 is a computation of mixing 
efficiency n , defined as the fraction of hydrogen, integrated 
over the flow cross-sectional area at a given station, that would 
react If complete reaction with the available oxygen were to occur. 
This parameter has been used for correlating cold-flow data (ref. 

9) and Is readily computed by COMOC as an output parameter by 
means of the Integration techniques utilized for measuring hydro- 
gen mass flow. Hence, figures 28 to 30 demonstrate that the 
virtual-source concept of transverse hydrogen Injection for the 
cold- flow configuration effectively simulates the Injection 
phenomenon . 

It Is therefore hypothesized that the virtual-source concept 
Is appropriate for combustion studies as well, and the experi- 
mental verification of this hypothesis Is sought. As a first 
step. It Is appropriate to measure the Influence of combustion 
on the virtual-source cold- flow configuration. Shown also In 
figure 30 are computations carried out to x^/d = 30 for the cold- 

flow simulation, where combustion of the hydrogen is allowed to 
occur according to the complete-reaction hypothesis (equation 
(35)). Note that the trajectory of maximum hydrogen mass fraction 
lies considerably above that for the cold-flow, non-reacting con- 
figuration. However, on the basis of mixing efficiency n , there 
Is very little difference In overall mixing between the cold re- 
acting and non-reacting cases. 

The cold-flow problem Is of marginal Interest, however, 
since the average equivalence ratio of the cold-flow configura- 
tion ((j) = 0.04) lies well below the design level for a practical 
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combustor. Note that equivalence ratio Is defined (ref. 11) 
strictly In terms of the global mass flows of hydrogen and air. 

For stoichiometric combustion, (j) = 1; for fuel-lean operation, 

4> Is less than 1. It Is coincidental that the cold-flow virtual- 
source configuration can be thermodynamically altered (only) to 
correspond to conditions similar to test point 4 of reference 10. 

To simulate the test configuration, the cold flow was computa- 
tionally vitiated by Imposing an arbitrary uniform background 
hydrogen concentration (of 1 percent) and augmenting the oxygen 
level of the base flow such that the corresponding composition 
of the computational test gas simulates the hot (wet) air used 
for the experiment. The total temperature of the computational 
simulation was approximately 2000 K; the corresponding mass flow 
of cold hydrogen for the vitiated virtual source yielded cp = 0.5 
for the simulation. Shown In figure 30 are the trajectory of 
maximum hydrogen mass fraction, the elevation trajectory, and 
the lateral spreading rate for the vitiated virtual-source simu- 
lation of the test configuration. Note that the elevation tra- 
jectory of the hot-flow configuration follows very closely the 
cold- flow data, a result possibly of the cold-wall (T^ = 0.5T^) 

boundary condition used for the computational simulation. It 
may also reflect the somev/hat lessened lateral spreading rate 
for the vitiated reacting case, as shown In figure 30. Mixing 
efficiency was not computed for this vitiated combustion case. 
However, equivalence ratio, as a node point parameter, can be 
computed at any point In the solution domain. At the far right 
of the mlxlng-ef flclency curve In figure 30, the experimentally 
determined range of equivalence ratio for test point 4 (ref. 10) 

Is compared with the computed values. The presented computational 
values are In qualitative agreement with the experimental extremums 
on the center plane at the duct exit. Analysis of the specific 
conditions studied In reference 10 were beyond the scope of the 
current Investigation. It Is suggested that future efforts with 
the COMOC code be directed In this area. 

H 2 Injection from Struts - Recent thinking on design of scramjet 

combustors Indicate that", depending upon flight Mach number, two 
fuel Injection modes will be required. For example, an experi- 
mental model of a strut Injector system (fig. 31) Is currently 
under construction at Langley Research Center for evaluation of 
combustor performance as a function of Injection mode. Design 
of this device was augmented by an earlier experimental program 
Intended to evaluate the essential character of the two distinctly 
different Injection modes proposed for this type of combustor 
(ref. 11). Schematics of the perpendicular and parallel Injection 
struts that are associated with current design technology are 
shown In figures 32 a and 32 b for cases 3-3 and 3-4 of table 3 . 

Also shown are the virtual-source simulations of the proposed 
Injection mode showing the location of the discrete Injectors as 
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well as the orientation of the virtual-source within the computa- 
tional domain. Figures 32c and 32d represent an exploded view 
of the virtual source simulation for each injection mode. The 
flux of enthalpy, mass and momentum was taken to be zero across 
each symmetry plane (see figures 32c and 32d) . At the wall bound- 
ary, the no slip boundary condition was imposed on the velocity 
components; heat flux and mass transfer were assumed to be zero. 
Other wall boundary conditions could have been employed for the 
enthalpy, e.g., specified wall temperature. However, in this 
series of computations, our prime interest was to compare the 
mixing between two different strut geometries since the wall tem- 
perature effects are of lower order. Obviously, this would not 
be the case if heat loads transmitted to the wall were of interest. 
The experlm.entally reported pressure distributions were used to 
obtain the required pressure gradient distribution. At this time 
the analysis is limited to constant area ducts (combustors), 
hence, area changes with distance downstream had to be neglected. 
However, the major effect of the area changes were implicitly 
considered by use of the experimental pressure field. 

Shown in figure 33 is the trajectory of the maximum hydrogen 
concentration as a function of injection mode for the experimental 
results reported in reference 11. Note Indeed that the perpendic- 
ular injection mode promotes much stronger mixing and hence pro- 
duces a combustion process that proceeds considerably more rapidly 
than that corresponding to the parallel injection mode. Also 
plotted in figure 33 are the pressure distributions used for the 
computations, as well as the computer distribution of equivalence 
ratio on the center plane at x^/d = 150 across one-half of a Jet. 

Note that for the parallel injection mode, the range of computed 
equivalence ratio is tv/lce that of the perpendicular case, in 
qualitative agreement with the data ranges from reference 11. 
Furthermore, the experimental evaluation of the differences in 
the flame shape for the two injection modes (ref. 11) with respect 
to apparent mixing rate is in agreement v;ith the maximum hydrogen 
trajectories by the virtual-source solution. 

Effects of Turbulence and Area Constraints - A number of investi- 
gators have shovm that it is often important to consider the 
interaction between the turbulent and reacting flow, e.g., ref. 

31. Mixing on the molecular scale is required for a chemical 
reaction to occur, hence, although the fuel and oxidizer may be 
mixed in the turbulent sense, l.e., on the macroscale, only a 
fraction of the fuel-oxldl zer is available for reaction. Chung, 
ref. 31, has proposed a simple model to consider this reaction 
limiting phenomena wherein the rate at which the reaction is 
allowed to occur is limited by the rate of turbulence dissipation. 
The effect of the turbulence-chemistry Interaction on the predicted 
flow field was analyzed using the virtual-source concept, the k+d 
turbulence model, and the reaction-rate limiting model of ref. 31. 
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Case 3-5 represents a case where the flow does not react, i.e., 
and 0 ^ mix on the microscopic scale and does not react. For 

this case the wall temperature was assumed constant at 296 . 0°K, 

hence, there Is a negligible amount of heat flux between the gas 
stream and wall. In cases 3-6 and 3-7 the H 2 /O 2 mixed on the 

molecular scale are assumed to autoignite. Table 4 shows the 

predicted surface heat flux at x^/d = 30 across the flow field 

for cases 3-6 and 3-7 where It is seen that Including the turbu- 
lence-chemistry Interaction In the model has reduced the peak 
flux by approximately 36 percent. Results of the nonreacting 
case ( 3 - 5 , flow conditions Identical to case 2-1) were used to 
examine the effect of the reaction on H 2 mixing. It was noted 

that the mixing efficiency n for case 3-5 was 76 percent, whereas 
when considering complete reaction (case 3-6) ri = 8l percent. 

This result suggests that for these conditions reaction does not 
have a significant effect on the mixing efficiency which Is In 
agreement with observations reported In reference 10. 

Three final cases (3-8a, 3-8b, and 3-9) were considered to 
examine the effect of Including the fact that the flow field Is 
ducted, and hence, generates a longitudinal pressure gradient 
due to boundary layer displacement effects, heat release, and 
friction losses. Rogers (ref. 10) performed his experimental 
studies In a 23 cm square test section divided by a splitter 
plate from which H 2 was Injected from a row of orifices. Since 
reacting flow was "^not of Interest In these (ref. 10) studies, 
the equivalence ratio was very low, i.e., approximately 0.04. 
Therefore, It Is not surprising that the primary result from these 
cases was that considering the ducting of the flow had no signi- 
ficant effect on the velocity, temperature, and H 2 concentrations 

for physically reasonable duct heights. Specifically, the pres- 
sure Increased by approximately 25 percent In cases 3~8a and 3-9 
with a slight drop In the maximum velocity (~ 2 percent). How- 
ever, the reacting flow could be forced to thermally choke by 
continually decreasing the duct height (Increasing mixture ratio 
to 0 . 50 ) to very small values, e.g., case 3-9 where X 2 ) ~ 2.0 cm. 

(compared to 13.5 cm which Is the half height of the wlnH tunnel 
section used In reference 10). 

The results presented for the cases of table 3 show that the 
analytical predictions are In general agreement with experimentally 
observed trends. Further use of the theoretical model Is definitely 
required. In particular In modeling data for other flow conditions, 
e.g., refs. 26, 27, 32. Areas of particular Interest which have 
not been considered In any significant detail are comparing model 
predictions with wall heat flux/temperature data, combustor pres- 
sure distributions, and combustor exit conditions. The utility 
of the COMOC code as a design tool will be significantly enhanced 
as the data base of flow conditions modeled using the code expands. 
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CONCLUDING REMARKS 


An analytical Investigation has been presented on the turbu- 
lent mixing and reaction of hydrogen jets Injected from multiple 
orifices transverse and parallel to a supersonic alrstream. The 
primary conclusions of this study are: 

(1) The analysis has Immediate utility In evaluating the 
mixing effectiveness of transverse Injection data since It has 

been tested In Its ability to model this type of data. 

(2) Turbulent mixing length theory, constant effective 
Prandtl number, and a Lewis number of unity provide reasonable 
agreement with transverse Injection data downstream of the 

near Injection region. Improvements In the correlation between 
predictions and experiment were obtained using a turbulent Prandtl 
number model and a two-equation turbulence model. 

(3) The effect of turbulence-chemistry Interaction could 
reduce the wall heat flux by as much as 36 percent. Comparisons 
with data are required to substantiate this predicted result. 

(4) Results demonstrate that the theoretical model embodied 
In the COMOC code has the potential for providing useful engi- 
neering design guidance for data covering a broad range of flow 
conditions, combustor geometries, and Injection models. 

(5) As with all complex turbulent flows, reliable theoretical 
predictions are only attainable upon testing a model's ability to 
accurately characterize the class of flows of Interest. Future 
efforts should be directed toward using the code In modeling data 
for a wide range of flow conditions. 

Although efforts herein were primarily directed toward 
comparing predicted and experimental mass fractions, other 

types of data are available, e.g., wall temperature and heat flux 
data, nozzle exit conditions, combustor pressure variations. The 
potential of developing COMOC as a design tool will be signifi- 
cantly enhanced as the data base of flow conditions modeled expands. 
In addition, areas for Improvement In the model will become appar- 
ent. Ultimately, the code may well both supplement and In some 
cases replace the costly experimental evaluation of prototype 
designs . 
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APPENDIX 


DATA DECK PREPARATION 


A description of how to prepare a data deck for the 3DBR 
variant of COMOC was reported In reference 5? In that version 
of the code, flows were considered where the pressure field was 
given and no area constraints Imposed on the flow. In modifying 
the code to consider constant area ducted flows ( 3DPNS-C0M0C ) 
some minor modifications to the data input procedure were required. 
Table A-1 Is a listing of the data deck used to run case 3-8b 
which Is a reacting, ducted flow. Other problems with different 
geometries or flow conditions can be treated using this data 
check as a starting point since approximately one-third of the 
COMOC data deck Is associated with standard call sequences as 
well as output format specifications and arrangement Instructions. 
If only scramjet combustors were to be analyzed then these 
standard data cards could be incorporated into the main program, 
and hence, removed from the data deck. They have been retained 
In order to allow for future modifications and additions for 
those Interested In expanding the analysis capability of the code. 
It Is suggested, however, that the standard data not be altered 
without consulting reference 7. Details on data deck preparation 
are presented below. The input may be conveniently divided Into 
four distinct data sets. Input data set 1 deals with the defini- 
tion of reference conditions, number of dependent variables, 
type of flow (reacting or non-reacting), etc. Input data set 2 
deals with the finite element discretization, wherein the number 
of nodes and finite element sizes and location are defined. 

Input data set 3 defines the type and format of the output, 
whereas data set 4 prescribes the initial conditions, pressure 
field and boundary conditions. The user of COMOC does not re- 
quire becoming familiar with input data set 3 unless the code Is 
modified to consider dependent variables other than u^ , U 2 , u^, 

H, , Yq , and Yj^ . In order to make such a modification, a 

familiarity with details of the codes internal structure Is 
required and the Programmer's manual should be consulted (reference 
7) . 


Data Set 1 

Cards 1 to ^ of Table A1 Represent Input Data Set 1 . 

Card 1.1 Starting In Column 1 FEEL used to start execution of 
COMOC 

Card 2.1 Starting In Column 1 Include 3DPNS If a constant area 
duct Is to be considered and hence pressure will be 
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Card 3.1-3. n 


NEQKNN 

IGAS 0 
1 

IFR 0 


1 

KDUMP 0 
1 

NPVSX 
NSCX 0 
1 

NSCY 0 
1 

NE1E2 0 
1 

Cards 4,l-4.n 


UINF 

T0PINF 

REEL 

T0 

TD 


computed internally. If the static pressure field 
is to be given then starting in Column 1 Include 
3DBR 

Namelist NAMEOl. A complete list of all NAMEOl 
variables and their default values is listed in 
subroutine FENAME. As seen from Table A1 of the 
one-hundred and sixty-seven parameters which can 
be read in NAMEOl, only eight need be considered 
for scramjet problems of Interest. The remaining 
parameters have been retained to facilitate future 
expansion and modification of the code. The seven 
Integer variables, which are required, are defined 
as follows. Note that in inputting in namelist 
form there is no designated column or order for 
each variable. The only requirement is that a 
comma follow Immediately after each integer value. 

Number of dependent variables to be Integrated in 
XI direction 

Isoenergetic flow with constant c 
General flows ^ 

Equilibrium composition or complete reaction 
(determined by version of subroutine GAS used 
in COMOC) 

Frozen composition 

Suppress debug output in gas subroutine 

Print debug output 

No. of entries in pressure table 

Uniform X3 interval in discretization 

Non-uniform X3 Interval in discretization 

Uniform X2 interval in discretization 

Non-uniform X2 Interval in discretization 

Laminar Flow 

Mixing length theory 

Namelist NAME02. A complete list of all NAME02 
variables and their default values is given in 
Subroutine FENAME. As seen from Table A1 of the 
one-hundred and seventy-four parameters which 
can be read in NAME02 only fifteen need be con- 
sidered and are defined below. The input mode 
is identical to NAMEOl. 

Reference (freestream) velocity (F/S) 

Reference stagnation temperature (°R) 

Reference length (P) 

Initial XI station (F) 

Length of XI solution, starting at T0(F) 
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DELP 

■VSTART 

XSCALE 

YSCALE 

CON 

XLAM 

PR 

SCT 

CVU 

Data Set 2 
Cards 5 to 
Card 5.1 

Card 6 . 1 
Cards 7.1-7 


Mode 1: 


Percent increment of TD at which output is desired 
Percent of TD at which transverse velocity (U2) 
computation starts - Default = 100%. 

Multiplier on Xo to convert discretization to feet 
as input under ^Mode 2 (see lines 7.1-7.n) 
Multiplier on X 2 to convert discretization to feet 
as input under Mode 2 (see lines 7.1-7.n) 

Constant in mixing length theory 
Constant in mixing length theory 
Effective Prandtl number 
Effective Schmidt number 

Equal to 0.009659 if simplified GPAHFT subroutine 
is employed, i.e., where simple reaction model 
is used. If equilibrium composition is employed 
use default value (do not input CVU). 


of Table A1 represent input data set 2 

Starting in column 1, FEDIMN calls the subroutine 
to generate vector lengths and array entry points. 
This card and its position cannot be altered by 
the user. 

Starting in column 1 LINKl - then after leaving 
a blank in column 6 - place a 1 in any column 
from 7 to 72. This card Instructs the program 
to call subroutine LINKl to read the finite 
element discretization In the (X2, X3) plane. 

n This set of cards can be input in one of two modes 
depending on the Integer values of NSCX and 
NSCY read in NAMEOl. Their function is to define 
the discretization. In either mode, the input 
card format is free - meaning no specific column 
has to be filled. A blank field indicates that 
a value has been read and the code search for the 
next value, which can be on the source card or 
the next card. A "T" represents the end of the 
data set. The two possible input modes are as 
follows . 

Automatic Uniform Discretization 
Occurs for NSCX = NSCY = 0 (in NAMEOl) where user 
set XSCALE = desired element width in the X3 
direction and YSCALE = desired element height in 
the X2 direction (in NAME02). On card 7.5 in any 
column indicate the number of the first node 
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(always 1). Skip one or more places and add the 
last node number In the X2 direction. Place a 
comma immediately after the last X2 node number. 
Repeat for X3 direction, e.g., if YSCALE = 0.004 
and X SCALE = 0.002 and we desired 21 x 2 nodes 
and 2 columns then card 7.5 would appear as 
1 21 , 1 2 , 

The sequence would be completed by placing a "T" 
either on this card or on a following card in any 
column from 1-72. The elements would be 0.004P 
high by 0.002P wide. 

Mode 2: Automatic Non-Uniform Discretization 

Occurs for NSCX = 1 and NSCY = 1. Set X3 dis- 
cretization first, X2 discretization second. 

Data are used in sets of 3 Integers at a time. 
First Integer identifies finite element Interval 
concerned, next two indicate element width (or 
height) as a dimensionless ratio, e.g., 3 1200 e 
3/1200 

e.g., 1 3 1200,2 1 600,3 5 1200,... 

T 

1 1 600,7 1 600,8 7 1200 , . . . 

T 

1 11, 1 4, 

T 

This generates a finite element discretization of 
11 node rows x 4 node columns. The element widths 
(Intervals between node columns) are respectively 
3/1200, 1/600, ... The height of the first 7 
element rows is uniformly 1/600, eighth is 7/1200, 
etc. These dimensionless Intervals are converted 
to feet by being multiplied by XSCALE and YSCALE 
(Cards 4. 1-4. 9) 

Data Set 3 


Cards 8 to 29 of Table A1 represent input data set 3 


Card 8 . 1 
Card 9.1 
Card 9 • 2 


Starting in column 1 COMTITLE, which designates 
that the following card will be a title card. 

Title card (columns l-80) which is printed on 
cover page of output. 

Starting in column 1 DONE indicating that this 
sequence is complete. 


Card 10.1 Starting in column 1 DESCRIPT - then in column 11 

or more - 204 followed by one or more spaces and 
a ”T" designating the end of this card. Note 
that any message can be placed after the "T" 
without affecting the data. 
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Cards 11.1- 
11. n 


Card 12.1 

Cards 13.1- 
27.1 


Cards 13.1- 
17.1 

Cards l8.1- 
19.21 


Cards 20.1- 
21.7 


Cards 22.1- 
23 .^ 


Up to ten title cards can be input which will be 
printed at each output station (XI) (hence n £ 10). 
These titles will head the generated output 
sequence . 

Starting in column 1 DONE 

These cards should not be changed without consulting 
the programmer’s manual (reference 7)- Only a 
very brief description of their function Is given 
below . 


Provide the output data heading descriptors. 


Designate what multiplying factor should be used 
to convert each output variable to units com- 
patible with the headings given by cards 13.1 to 
13.12. 

Designate vector location of variables to be 
printed per headings given by cards 13.1 to 13.12 
after being multiplied by factors designated by 
cards l8.1 to 19.21. 

Designates what variables are to be printed out 
at all node point locations . Any variable ending 
with the numerals 248 represents an dependent 
variable. The numerals preceding the 248 designate 
what dependent variable is to be printed. For 
Table A1 conditions the following number-parameter 
relations hold. The variable location contains 
the dimensionless value and Is printed as such. 

Free format Is used. 


Location 

Variable 


1248 

^1 


285 

static temperature 


320 

static enthalpy 


284 

density 


10248 

N 2 elemental mass 

fraction 

3248 

^3 


278 



1248 

stagnation enthalpy 

9248 

H 2 elemental mass 

fraction 

8248 

0^ elemental mass 

fraction 


I 
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1247 

^eff 

334 

^^^eff 

292 

y (laminar) 

314 

^'^^eff 


Cards 24.1- 
25.2 

Card 26.1 

Cards 28.1- 
28.4 

Card 29.1 
Data Set 4 
Card 30.1 


Card 31.1 


Card 32.1 
Card 32.2 
Card 33.1 


Indicates the multiplier to be used on the variables 
designated by cards 22.1 to 23.4 prior to print 
out. Here all fourteen variables are scaled by 
unity, l.e., they retain their dimensionless form. 

Starting In column 1 COMOC to key code that the 
problem description Is to be provided on the 
following cards. 

Problem description on four cards columns l-80. 

This description Is printed once at start of run. 

Starting In column 1 DONE 


Starting In column 1 ’VX3ST designating that the 
next card will give the XI locations at which 
pressure will be Input - XI In feet and pressure 
In PSFA. If 3DBR option Is used this pressure 
profile will be used for the complete calculation. 
If 3DPNS Is used this pressure profile Is used 
to specify the Initial pressure and to compute 
the Initial pressure gradient. 

Starting in any column, the XI locations for 
pressure data. Each value must be separated by 
one or more spaces. No Integer value Indicating 
the number of points to be read Is required. The 
code counts the values on the card. Data set is 
complete by placing a "T" after the last value. 

Any number of cards may be used. 

Starting In column 1 VPVSX 

Pressure points In PSFA. Same format as card 31.1. 

Starting In column 1 IPINT and In any column after 
10 Include -1. This card designates that the 
Integer array numbers of the dependent variables 
will follow on card 34.1. The program will 
integrate the first NEQKNN of these values and 
also U 2 . 
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Card 3^.1 


Starting '±n any column list the integer values 
of the dependent variables which correspond to 
the following key: 


Card 35.1 


Mode 1: 


Integer 

1 

2 

3 

5 

6 

7 

8 

9 

10 


Variable 


u 

u 


1 

2 


u 


3 


H 


Open 

Open 

Open 



2 



A free format is used, hence, the sequence is 
terminated with a "T." The open Integers 5, 6, 
and 7 have been retained to accommodate the 
Introduction of additional variables such as 
turbulence kinetic energy. 


Boundary conditions are required at each boundary 
node point for each dependent variable. The 
following input instructions refer to any 
dependent variable. The data cards defining 
boundary conditions for other variables follow 
the preceding boundary condition data cards, 
l.e., cards 35.1 and 36.1 are repeated for each 
dependent variable. 


There are three modes under which boundary condi- 
tions at each boundary node may be defined: 

(1) zero normal gradient, (2) value fixed to Its 
Initial value, and (3) normal gradient defined 
as a finite value, Eq. A1 . 


9n 


a 


(Al) 


Zero Gradient Boundary Condition 
This represents the default boundary condition 
hence any node whose boundary condition is not 
defined will automatically be given a zero 
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Mode 2: 


Mode 3" 


Card 36.1 
Mode 1: 
Mode 2: 


Mode 3 : 


2 

gradient boundary condition, e.g., — g— - = 0 Is 

generally valid on all boundaries, hence cards 
35.1 and 36.1 are not Input for this variable. 

Boundary Node Fixed to Its Initial Value 
Starting in column 1 KBNO and after column 10 
the number of the dependent variable whose 
boundary condition will be defined on card 36 . 1 . 

Boundary Node with Defined Normal Gradient, Eq. A1 
Starting In column 1 KBNO, in column 11 the 
number of the dependent variable, and In column 
21 the integer 1 which indicates that mode 3 
boundary condition will be imposed. 


Card 36.1 not used for this mode 

At each station the cross-sectional computational 
domain corresponds to a rectangular shape. Herein 
we define, for future reference, the boundaries 
as TOP, BOTTOM, LEFT, and RIGHT. The code has 
been programmed so that nodes along a given boundary 
are not required to have the same boundary condi- 
tion. 

To fix values along the boundary, the general card 
image is as follows: 


Column = 

1 11 

• • 

21 

31 

41 

1 — 1 1 

ir\\ • 

61 

11 

Card = 

• • 

• • 

BOTTOM NB 

TOP 

NT 

RIGHT 

NR 

LEFT 

NL 

This card 

designates 

that 

the 

first 

NB, 

BOTTCM 

9 


NT, TOP, NR, RIGHT, and NR, LEFT boundary nodes 
will be held constant. By leaving any of the 
Integers (NB, NT, NR, or NL) out all nodes are 
held constant. If only TOP and LEFT nodes are to 
be held constant, the card image would contain 
the word TOP in column 1 and LEFT in column 21. 

For those nodes not defined the default boundary 
condition (mode 1) applies. 

As in mode 2 we make use of the TOP, BOTTOM, LEFT 
and RIGHT designation of the boundaries. However, 
unlike mode 2 all nodes along a given boundary 
must be treated the same. Starting in column 1 
the descriptor of the Boundary is given followed 
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Card 37.1 


Card 38.1 


Card 39.1 


Card ilO.l 


Card ill.l 


Cards 42.1- 
42. n 


by the value of A starting in column 21 and 
followed by a blank space and then the Integer 2, 
e.g. , if A = . 342 
Column =1. 21 

Card = BOTTOM .342 2 

Similarly, conditions for TOP, LEFT, and RIGHT 
may be employed. The 2 following the value 
Instructs the code that the input boundary condi- 
tion will be input in dimensionless form, e.g., 
if f H u]_ then A must be defined consistent with 
the velocity non-dimenslonallzed by u^ ^ and the 
length scale normalized by 

Starting in column 1 PRINT which designates that 
the message appearing on the following card 
will be printed before execution of the run. 

Anywhere in columns I -80 a message which is 
printed before execution which has generally 
been used to describe the type of boundary 
conditions employed. 

Starting in column 1 LINK3 then in any column 
after 11 the Integer 4, which Instructs the 
program to call subroutine LINK3. The Integer 
4 Instructs the code to execute the DIMEN section 
of LINK3 (see reference 7). 

Starting in column 1 LINKl then in any column 
after 11 the Integer 3 which Instructs the 
program to call subroutine LINKl. The Integer 
3 Instructs the code to execute the GEOMEL 
section of LINKl (see reference 7). 

Starting in column 1 VTBVIP - then in any column 
after 11 - the Integer -58. VTEMP designates 
that the stagnation temperatures at each node 
point will be input on the following cards. The 
Integer -58 Instructs t^e code to put the 
temperatures (Input In ^R) into dimensionless 
form by dividing by variable number 58 . 

Starting In any column the stagnation temperatures 
at each node point. The order that the values 
appear on the card(s) correspond to the node 
number. The nodes are numbered as follows. 
Consider a problem with n rows and m columns. 

Node 1 corresponds to the lower left corner node 
and nodes 1 to m designate the nodes along the 
bottom. The second row from the bottom is 



numbered from m+1 to 2m and so on with the upper 
right node number being nxm. If N successive 
values X are identical then these values may be 
input as N*X rather than X space X space . . . X . 

The end of data is designated by a "T" after the 
last entry. 

Card 43.1 Cards 43 » 44, and 45 are used to input one of the 

velocity components in ft/s. These three cards 
must be repeated for each of the three velocity 
components u^, Ug, and u^. If a velocity is not 

defined, it is set equal to zero. 

Starting in column 1 VYY then in any column after 
11 the Integer -27. VYY designates that one of 
the dependent variables to be Integrated forward 
will be input on the following cards. The 
integer -27 instructs the code to divide the input 
variable by the value stored in RARRAY(27) which 
is the reference velocity. 

Cards 44.1- Value of the U]_ velocities in ft/sec using same 
44. n free format mode as used in cards 42. 1-42. n. 

Card 45.1 Starting in column 1 VYYEND - then in any column 

after 11 - the Integer 1. This card indicates 
that the variable input above (cards 44. 1-44. n) 
correspond to values of dependent variable 1 (u^). 
Repeat cards 43-45 for U 2 and u^ and designate 

what velocity component is being input by the 
Integer on card 45.1. 

Card 46.1 Starting in column 1 VYY. Note the absence of any 

Integer Indicates that the variables are in the 
proper dimensionless form. 

Cards 47.1- Values of dependent variable at each node point. 

47. n 

Card 48.1* Starting in column 1 VYYEND then in any column 

after 11 the Integer 9- This card Indicates 
that the variable input above (cards 47.1-47.n)„ 

“2 

correspond to values of dependent variable 9 (Y ). 

Cards 49.1- These cards are not to be altered. They direct 
56.1 the code to the integration loop and control the 

printout of reference values which appear at the 
end of the output description. 

*Inltlal total enthalpy is obtained from stagnation temperature. 

If initial values of O 2 and N 2 are not input, th.ey are computed 
Internally using the assumption that the mass other than H 2 is 
23.2% O 2 and 76 . 8 ^ N 2 j l.e., air. 

47 


I 



J=- 

00 


TABLE A1 

COMOC SAMPLE INPUT DATA DECK 


Data Set 1 (1.1 to 4.6) 


1 — i 
( — 1 

FEBL 





2.1 

3DPNS 





3.1 

FENAME 





3.2 

CNAMEOl 





3.3 


NPVSX=2, 

NSCX=i , 

NSCY=1, 


3.4 

NE1E2 

= 1 , 




3.5 


NEQKNN=3, 

IGAS=1 , 

IFR=1, 

KDUMP=l 

3 . 6 

£,ENO 





4.1 

&NAMF02 





4.2 


UINF=2272. , 

T0FINF=533.0, 

REFL=. 003333333, 


4.3 


XSCALE=0. 003333333, 

YSCALE=0. 00333333 3, 

VSTART=101.0, 


4.4 


T0=0,0, 

TD=0. 1, 

DELP=05 .0, 



4.5 CnN=0-435, XLAM=0.07, PR=0.7, SCT=G.7f CVU=C. 009659 , 

4.6 -SEND 

Data Set 2 (5.1 to 7.6) 

5.1 FEOIMN 

6.1 L INK 1 1 SETUP 

7.1 1 75 100, 2 50 100, 3 125 100, A 150 100, 5 225 100, 

7^2 T INCREMENTS BETV«EEN X3 , NODE-NUME R A TCR-DE NOM I N ^T0R 

7.3 1 1 4, 2 1 4, 3 1 2, 8 1 2, 9 1 1, 10 9 2, 11 27 2, 12 55 2, 

7.4 T INCREMENTS BETWEEN X2 

7.5 I 13, 1 6, 

7.6 T 13 ROWS AND 6 COLUMNS NORMALIZED BY LREF, HENCE X-Y SCALES =LREF 


TABLE A1 (Contd) 


Data Set 3 (8.1 to 29.1) 


8.1 COMT ITLt 

9.1 VIRTUAL SOURCE - 30PNS - MLT - CCPf-LETE REACTION 

9.2 DONE 

_,0.1 DESCP IPT 204 T DESCRIPTIVE TITLE AT REGIMNG OF OUTPUT. 

11.1 VIRTUAL SOURCE - 30PNS - MLT - CCNFLETE REACTION 

11.2 

^2.1 done 

13.1 DESCRIPT 332 T lOPAR PARAMETER TITLES FCP OUTPUT, 

13.2 REFERENCE ENGLISH-FT ENGLISH-IN R-K-S C-G-S 

x3.^ LENGTH FT .IN M CM.... 

13.^ VELOCITY FT/S N.A... M/S CM/S..., ... 

13.5 DENSITY L8M/FT3 N.A KG/ M3 G/CC 

x3.6 TEMPERATURE RANKINE.... .N.A., KELVIN..... .N.A 

13.7 ENTHALPY BTU/LBM,... ,n,a KJ/KG. .N.A........ 

13.8 FROZ.SPEC. HEAT .BTU/LBM-R,. .N.A KJ/KG-K N.A 

s. 3.9 VISCOSITY...... .LBM/FT-S... .N.A NT-S/M2 POISE 

13.10 LOCAL PRESSURE .PSF .PSI........ .NT/M2 .TQRR 

13.11 local solution MACH NO. DPDX 1 ( LRF/ FT 3 ) MAX. H2 CUNC. MIX EFF.(ETA) 

13.12 Xl/LREF OXl/LREF EFSILCN DXIMIN/LREF 

14.1 DONE 

15.1 OESCRIPT 203 T IFMTHD TITLES FCR OUTPUT DEPENDENT VARIABLES. 

16.1 Ul/UREF T/TREF HSTAT/FPEF RFC/RHORFF ELEM.N2 MAS.FRAC 

16.2 U3/UREF CPF/CPFREF HTQT/HREF ELEM.H2 MAS.FRACELEM,G2 MAS.FRAC 

16.3 eff.mu/muref eff. prandtl no.mu/muref eff.schmict no. TKE/EKNINF 

diss/epsinf 

17 1 done 


VT) 



o 


TABLE A1 (Contd) 


18.1 MPARA -I 

19.1 5*2 

19.2 2 2 162 164 163 

19.3 222 164 163 

19.4 222 170 174 

19.5 22 2 165 2 

19.6 2 -175 2 2 2 

19.7 2 2 2 176 2 

19.8 22 2 177 178 

19.9 2 2 169 168 167 

19.10 2 2 2 2 2 

10.20 2 2 2 2 

19.21 

^0.1 fONUMB -1 

21.1 gqq 

21.2 5*200 999 

21.3 200 4*43 200 27 200 2*27 200 1C 200 2*1C 200 5B 200 58 200 

21.4 200 97 200 97 200 200 30 2C0 30 200 200 38 200 2*38 

21.5 999 

21.6 39 4*36 200 61 100 134 122 11 12 14 85 

21.7 

22.1I0SAVE -1 

23.1 1248 285 320 234 10248 

23.2 3248 278 4248 9243 8248 

23.3 1247 334 292 314 

23.4 T UfT,HS,RHn»N2,VfCP,HT0T,H2,n2,DIFU,PR NC . , LA . V I SC . t S C T . NO . 

24.110MULT -1 

25.1 14*2 

25.2 T U,T, HS,RHQ,r42,V,CP,HTQT,H2,O2,0IFU,PR NC. ,LAM. VI SC. , SCT .NO . 

26.1COMOC 

27 .IDFSCRIPT 

28.1 VIRTUAL SOURCE - 3DPNS - MLT - CCN'PLETE REACTION 

28.2 TURBULENCE MODEL EMPLOYED IS DESCRIBED IN USER S MANUAL NASA CR- 1 3245 0 , 1974 . 

28.3 CALCULATIONS ARE STARTED USING VIRTUAL SOURCE CCNCEPT TC REPLACE 

28.4 COMPLEX NEAR INJECTION FLOW FIELD. 

29. InONE 



TABLE A1 (Contd) 


Data Set 4 (30.1 to 56.1) 


3U.1VX3ST 
31.1 0.0 

100. 

T XI TABLE 

FOR PRESSURE 

32.1VPVSX 
32.1 193. 

193. 

T PRESSURE 

TABLE PSF 

33.1IPINT 

-1 




34.1 1 A 9 8 10 2 3 T 

35.1K8NO 1 

36.1BOTTQM DONE 

37.1PR INT 

38.1 FIXES LI (VARIABLE NC. 1) ALONG WALL TO INITIAL VALUE 

39.1LINK3 4 

40.1LINKI 3 

41.1VTEMP -58 

42.1 78*533,0 T 

43.1VYY -27 

44.1 6*0.0 

44.2 6*795.0 

44.3 6*1503. 

44.4 6*1660. 6*1759. 

44.5 2*1550. 4*1833. 

44.6 3*1550. 3*1892. 

44.7 2*1550.0 4*1942. 

44.8 1550. C 2*2272.0 3*1985.0 

44.9 3*2272.0 3*2068.0 12*2272.0 

44.10 6*2272.0 

44.11T INITIAL U1 PROFILE 
45.1VYYENO 1 


DIMEN 

GEOMFL 


VJl 



VJl 

TABLE A1 (Contd) 


46 . IVYY 

47.1 30+0.0 2+1.0 A + C.O 3+i.O 3*0.0 2+1.0 4 + 0.0 1 

47.2 6+0.0 

47.3 T INITIAL H 2 MASS FRACTION PROFILE 
48 . 1 VYYEND 9 

49.IQKMNT 
5O.IDESCR IPT 
5I.IDONE 
52.IDESCRIPT 3 

53.1 REFERENCE LENGTH, LREF 

53.2 REFERENCE VI SCO S I TY.L AM I NAR VALUE 

53.3 EVALUATED AT REF. TEMPERATURE. 

53.4 FREESTREAM VELOCITY AT XO(=UREFJ 


53.5 STAGNATION TEMPERATURE ( CON S TANT , = TREF ) 

53.6 FREESTREAM DENSITY AT XO(=RHCREFJ 

53.7 FREESTREAM MACH NUMBER AT XO 

53.8 STATIC PRESSURE AT XO 

53.9 NUMBER OF NODES 
53.IONUMBER OF FINITE ELEMENTS 
54.IDONE 

55 . 1 END 
56.IEX IT 


0 23+0.0 


4 2 FT. 

36 LBM/FT-S 
27 FT/S 


56 DEG R 
1C LBM/FT3 
154 
9 PSF 
-16 
-14 
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Table 1 

Sensitivity of Predictions to Various Parameters 


A. Conditions Used in all Cases 

^r 

d. In 

s/d ..... 

Distance from H2 Injection Port to Opposite Wall 

Free Stream Mach Number 

Stagnation Temperature (Primary), °K 

Stagnation Temperature 

Equivalence Ratio (m„ /ih . )/.0291 

rl 2 0,1 r 

Type of Plow 

^Simulates Injection Into an Infinite Air Stream 


1.0 

.04 

12.5 

13-5 cm 

4.03 

300° 
300° 
< . 01 * 

Non Reacting 


B. Purpose of Run 

Case Purpose 

1-1 Evaluate mixing length theory 

1-2 Double Discretization 

0 

1-3 Investigate effect of variations in turbulence constants X and N^^ 

1-4 Evaluate effect of finite cross flow (u^) velocity 

1-5 Investigate effect of turbulent shear stress model using tensorlal 

eddy viscosity 

1-6 Evaluate two equation turbulence model 

1-7 Evaluate turbulent Prandtl number model of ref. 30 



Table 2 


Data Used to Study Ability of Mixing Length Theory to Model 

Injection Data 

A. Conditions Used in all Cases 


Effective Prandtl Number O. 9 O 

Equivalence Ratio <0.01 

Free Stream Mach Number 4.03 

Stagnation Temperature 

(H^ and Air) 300°K 

Type of Plow Non Reacting 

B. Cases Considered* 


Case 

q.r 

s/d 

Ref. 

2-1 

1 

12.5 

9 

2-2 

0.5 

12.5 

9 

on 

1 

CM 

1.5 

12.5 

9 

2-4 

1.5 

6.25 

9 

2-5 

0.5 

6.25 

9 

2-6 

1.0 

00 

8 

1 

CM 

0.5 

CO 

8 

2-8 

1.5 

CO 

8 

*Note : 

Case = 1.0, s/d = 

6.25 not considered 

due to 


absence of symmetry In data. 
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CO 


Table 3a 


Cases Used to Study (a) Virtual Source Concept (b) Reacting Plow (c) Ducted Plow"'" 


Case 

^r 

s/d 

h/d 

M 

CO 

<P 

Mixing 

Model 

Reaction 

Model* 

3-1 

1.0 

12.5 

oo 

4.03 

O 

O 

4:=- 

MLT 

II 

3-2 

1.0 

12.5 

00 

4.03 

0.04 

MLT 

(air) 

II 

on 

1 

on 

- 

8.93 

10.0 

2.25 

0.60 

MLT 

(vitiated) 

II 

3-4** 

1.2 

8.92 

8.0 

2.25 

0.62 

MLT 

II 

3-5 

1.0 

12.5 

OO 

4.03 

0.04 

k+d 

I 

3-6 

1.0 

12.5 

00 

on 

o 

O 

o 

k+d 

II 

3-7 

1.0 

12.5 

00 

4.03 

0.04 

k+d 

III 

3-8a 

1.0 

12.5 

100 

4.03 

0,04 

MLT 

I 

3- 8b 

1.0 

12.5 

100 

4.03 

0.04 

MLT 

II 

3-9 

1.0 

12.5 

18 

4.03 

o 

o 

MLT 

II 


All cases used virtual source concept, cases 3-1 and 3-5 are non-reacting, cases 3-8 
and 3-9 are ducted, all others have a free boundary. 

*See Table 3b 

**Case 3-3 Is parallel Injection from strut, case 3-^ is perpendicular Injection; d 

Is Injector throat diameter (0.15 In. for 3-3, 0.187 for 3-^) otherwise d Is- sonic 
orifice diameter. 



Table 3b 


Reaction Models* 


I. 

Frozen 

Flow : 

(a) 

^2 

+ 

(b) 

°2 


(a) H 2 + (b) O 2 


II. A. 

(b) > 

(2a): 

(a) 

^2 

+ 

(b) 

°2 


(a) H 2 O + (b - |a) O 2 


II. B. 

(b) < 

(2a): 

(a) 

H2 

+ 

(b) 

°2 


(2b) H 2 O + (a - 2b) H 2 


III. A. 

(b) > 

(2a): 

(a) 

H2 

+ 

(b) 

°2 

-y 

(aR) H 2 O + (b - 1 aR) O 2 + 

(a(l-R)) H 

III.B. 

(b) < 

(2a): 

(a) 

^2 

+ 

(b) 

°2 


(2bR) H 2 O + (a - 2bR) H 2 + 

(b(l-R)) 0 


Notes: 1) ( ) represents number of moles 

2) R Is fraction of moles mixed on microscale and available for reaction 


ui 


Table 4 


X3/d 

0 

0.75 

1.25 
2.5 
4.0 

6.25 


Heat Flux Distribution at x^/d = 30 
for Virtual-Source Simulation* 


Case 3-6 
(Fully Reacted) 

V«w>r 
0.785 
0.853 
1.0 
0.230 
0.095 
- 0.015 


Case 3-7 
(Rate Limited) 

0.570 

0.610 

0.636 

0.108 

0.027 

- 0.014 


1) Fixed wall temperature, = 296. 0°K 

2) q, )^ = 0.33 Mw/m^ 

W I 


*Note : 



FOREBODY (INLET) 

(la) 

^ AFTERBODY (EXHAUST NOZZLE) 

HYPERSONIC PROPULSION SYSTEM SCHEMATIC COMBUSTOR LENGTH WITH 
SCRAMJET OPERATION ABOVE MACH 3 IN-STREAM FUEL INJECTION 



Figure 1 . NASA Hypersonic Vehicle 



(b) Top view 


Figure 2(a). Photographs of the Mixing-Reacting Flow Field of the Perpendicular-Injection Strut 






(b) Top view 


o^ 

U) 


Figure 2(b). Photographs of the Mixing-Reacting Flow Field of the Parallel Injection Strut 
















Boundary-Layer Thickness 5/6^ 
Mach Number Me 
Pressure 10"® 



Figure 4. Computed Supersonic Boundary Layer Parameters, 

M = 5, Rex = 0.83(5)/m, 13 = 0.5 65 
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Skin Friction 1/2 Cf x 10' 


















(a) Fully developed 


(b ) Developing 


Figure 6. Channel Flow Solutions Computed with the Parabolic Navier-Stokes Variant Equations 















Symbols are Best Symmetry Plane Fit for Data of Rogers (Ref. 5 ) 


Traverse Displacement, X2/d 



Lateral Displacement, X3/d 


Figure 9. Cubic Spline Interpolated Hydrogen Mass Fraction Contours 
for Single-Jet, qj- = 1.0, x/D =30 


70 











-<I 

J=- 


qf = 1.0; s/d = 12.5; Exp ref . 9; Predicted (MLT) 

Predicted (MLT, Twice Discretization) 



0 2 4 0 2 4 


Hydrogen Mass Fraction, Percent 


Figure 13. Comparison between Experimental and Predicted H 2 Mass Fraction Profiles along Center Plane (x 3 /d = 0) 
and Wall (x 2 /d = 0). Showing Effect of Doubling Discretization, Case 1-2 




<lr = 12.5 ; Exp ref. 9 ; Predicted (MLT, X= 0.09, Np^= 0.9) 

Predicted (MLT, X = 0.07, Np^ = 0.7) 



Figure 14. Comparison between Experimental and Predicted H 2 Mass Fraction Profiles along Center Plane (x 3 /d = 0) 
and Wall (X 2 /d = 0). Showing Sensitivity of Predictions to Turbulence Model Constants, Case 1-3 




Figure 1 5. Initial Cross-Stream Velocity Distribution (xj /d = 30) 
and Computed Distributions at Xj/d = 60 and 120 





= 1.0; s/d = 12.5; Exp ref. 9; Predicted (MLT) 

Predicted (MLT, U 3 =h 0) 



Hydrogen Mass Fraction, Percent 



0 2 4 0 2 4 


Hydrogen Mass Fraction, Percent 


Figure 16. Comparison between Experimental and Predicted H 2 Mass Fraction Profiles along Center Plane (x 3 /d - 0) 
and Wall (x 2 /d = 0). Showing Effect of Including Transverse Velocity, Uo, Case 14 




Qf = 1.0; s/d = 12.5; Exp ref. 9; Predicted (MLT) 

Predicted ( 2 ^ 3 ) 




Hydrogen Mass Fraction, Y^, Percent 
(a) ei 2 = e,e ,3 = 2e 




Hydrogen Mass Fraction, Y*^, Percent 
(b) e,2 = 2e, ei3 = e 


Figure 17. Comparison between Experimental and Predicted H 2 Mass Fraction Profiles along Center Plane 
(x 3 /d = 0) and Wall (x 2 /d = 0). Showing Effect of Tensorial Eddy Viscosity, Case 1-5 





qr = 1.0; s/d = 12.5; Exp, ref. 9; Predicted (MLT) 

Predicted (k + d) 



Hydrogen Mass Fraction, Y^, Percent 


Figure 18. Comparison between Experimental and Predicted H 2 Mass Fraction Profiles along Center Plane 
(x 3 /d = 0) and Wall (x 2 /d = 0), using 2 Equation Turbulence Model, Case 1-6 


VO 
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CO 
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= 1.0 ; s/d = 12.5 ; Exp ref. 9 ; — — — — Predicted (MLT, Np^ = 0.9) 

- Predicted (MLT, Np^. Model from Ref. 27) 



0 2 4 0 2 4 

Hydrogen Mass Fraction, Y*^, Percent 



Hydrogen Mass Fraction, Percent 


Figure 19. Comparison between Experimental and Predicted H 2 Mass Fraction Profiles along Center 
Plane (x 3 /d = 0) and Wall (x 2 /d = 0). Case 1-7: MLT, Npj. Model from Ref. 27 











Qr = 1.5 ; s/d = 12.5; 


Exp ref. 9; 


Predicted (MLT) 





Figure 21. Comparison between Experimental and Predicted H 2 Mass Fraction Profiles along 
Center Plane (x 3 /d = 0) and Wall (x 2 /d = 0). Case 2-3, MLT, N|j. = 0.9 






Hydrogen Mass Fraction, y'^, Percent 


Figure 22. 


Comparison between Experimental and Predicted H 2 Mass Fraction Profiles along 
Center Plane (x 3 /d = 0) and Wall (x 2 /d = 0). Case 2-4, MET, N|j- = 0.9 
















qr = 1.0; s/d = °° ; Exp ref. 8 ; Predicted (MLT) 



Hydrogen Mass Fraction, Percent 



Hydrogen Mass Fraction, Y*^, Percent 


Figure 24. Comparison between Experimental and Predicted H 2 Mass Fraction Profiles along 
Center Plane (x 3 /d = 0) and Wall (x 2 /d = 0). Case 2-6, MLT, Npj. = 0.9 






qj. =0.5 ; s/d = °° ; Exp ref. 8 ; Predicted (MLT) 
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Figure 25. Comparison between Experimental and Predicted H 2 Mass Fraction Profiles along 
Center Plane (x 3 /d = 0) and Wall (x 2 /d = 0). Case 2-7, MLT, Npj. = 0.9 
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Figure 27. Transverse Injection into a Turbulent Boundary Layer 
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Figure 28. Comparisons between Experimental and Predicted H 2 Mass Fraction Profiles 
at xj/d = 30 and Various X 3 /d Stations. Virtual-Source Concepts used 
to Start Calculations at xj/d = 0 
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Figure 29. Comparisons between Experimental and Predicted Mass Fraction Profiles 
at xj/d = 120 and Various X 3 /d Stations. Virtual-Source Concept used 
to Start Calculations at X I /d = 0 
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Figure 30. Transverse Cold Hydrogen Injection with Virtual-Source Simulation (qj- = 1 .0; s/d = 1 2.5 of Reference 9). 
Vitiated Reacting Flow Data of Reference 1 0 does not Exactly Correspond to these Conditions since 

Qi- = 1.26, s/d = 10.5, 0 = 0.63 









Parallel injection strut 





Maximum hydrogen mass fraction, percent 


A Perpendicular strut 

O Parallel strut 

Solid symbols denote data 




Figure 33. Analytical Evaluation of Two Supersonic Strut Injectors from 
Virtual-Source Simulation, tj = 3; <p = 0.6; s/d = 9 


